THE BELL SYSTEM 

TECHNICAL JOURNAL 

volume xlvi September 1967 number 7 

Copyright © 1967, American Telephone and Telegraph Company 

Statistical Analysis and Modeling of the 

High-Energy Proton Data From the 

Telstar® 1 Satellite 

By J. D. GABBE, M. B. WILK and W. L. BROWN 

(Manuscript received October 4, 1966) 

This paper deals with the analysis of data from the omnidirectional 
high-energy proton detector on the Telstar 1 satellite. The main accom- 
plishment is the development of relatively simple (empiiical) mathematical 
models which give a statistically accurate representation of the measured 
spatial distribution of intensity of protons with energies between 50 and 
ISO MeV. 

These models depend upon the fitting of 8 (or 9 or 10) coefficients based 
on samples containing approximately 1000 of the nearly 80,000 experi- 
mental observations. The nature of the model for the average omnidirec- 
tional counting rate permits its closed form transformation to the equivalent 
equatorial pitch angle distribution. 

Sufficiently accurate fits were achieved so that the residuals {equal to 
observed minus fitted) coidd be productively examined for possible depend- 
ence on variables other than the two magnetic coordinates used in the 
fitting. One consequence of this icas the detection of instrumental suscep- 
tibility to temperature and bias voltage changes, which led to an objective 
partitioning of the data. 

The present paper has several evolutionary aspects: In particular, a 
series of one-dimensional fits was employed as a base for developing a 
two-dimensional model; a preliminary analysis of all the data ivas used 
to guide the rejection of outliers; a first two-dimensional fit to all the data 
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led to a data-independent basis for partitioning the data; the mode oj 
selection of a sample oj data, to which the two-dimensional model was 
fitted, changed as deeper insight into the importance oj this issue developed; 
and, ajter a very satisfactory fit to the data was attained, the model was 
improved by specialization and reparameterizalion so as to overcome some 
statistical dejects and to achieve greater physical meaning. 

The data cover the time period between July 1962 and February 1963, 
and the spatial region bounded by 1.09 R.^R^ 1.95 R,andO^\< 58°. 
Flux maps having a relative accuracy oj about two percent are derived 
jrom the fit and presented. The temporal behavior oj the intensity is ex- 
amined and some changes are noted. The maximum value oj the omni- 
directional flux oj protons with energies between 50 and 130 MeV is jound 
to be [5.711:1] X 10 3 protons/cm 2 sec at L = 146 on the magnetic 
equator, in good agreement with other experiments. Relative flux values 
and energy spectra are consistent with the generally accepted picture oj 
the proton distribution. 
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I. INTRODUCTION 

This paper deals with the analysis of data from the omnidirectional 
high-energy proton detector on the Telstar® 1 satellite. The main ac- 
complishment is the development of a relatively simple (empirical) 
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mathematical model which gives a statistically accurate representation 
of the measured spatial distribution of protons with energies between 

50 and 130 MeV. 

The Telstm-® 1 satellite was launched into a 45°-inclined orbit with 
an apogee of 5600 km and a perigee of 950 km on day 191 (July 10) , 

1962. The period of precession of the apsis was 180 days. The satellite 
was instrumented to measure fluxes of energetic particles; in particu- 
lar, counting rates of protons with energies above 50 MeV were re- 
corded. Two thousand hours of telemetry was received during the ac- 
tive life of the satellite, which terminated on day 52 (February 21), 

1963. The satellite and associated systems have been described in de- 
tail. 1 The particle-detection instruments have been documented 2 and 
some of the experimental results have been presented. 3, 4 - B 

The above-mentioned presentations of information concerning the 
earth's radiation belts have been principally graphical in format, ow- 
ing to the complexity of the belts and the limited understanding of the 
details of the processes affecting them. 

An accurate analytical representation of the data would enable con- 
venient interpolation, extrapolation, and transformation. Thence it 
would be practical to make extensive comparisons with the results of 
other experiments and with various theoretical predictions and to sum- 
marize, analytically, such features as the equatorial omnidirectional 
counting rate and the approximate size of the equatorial loss cone. In 
addition, an empirical mathematical model would facilitate the study 
of temporal fluctuations in various regions of space. Of course, a good 
analytical representation, even though empirical, may also stimulate 
deeper physical insight and theories. 

The present study was directed toward the development of a math- 
ematical function which would, when fitted to the data, provide a con- 
venient, concise and precise summary description. The mathematical 
model (s), which are herein presented, were empirically evolved, using 
the knowledge that the intensity distribution of these protons is, in 
the main, not rapidly variable in time. Even more specifically, the 
assumption has been that fluctuations in observed counting rates at a 
fixed point in space relative to the earth are independent random vari- 
ables. Further, the main effort of the analysis has been to try to relate 
the observed counting rates to a two-dimensional magnetic coordinate 
system derived from three-dimensional spatial coordinates by mapping 
the known earth's magnetic field onto the field of a magnetic dipole. 

The mathematical models which are used depend upon fitting of be- 
tween 8 and 10 coefficients based on samples containing approximately 
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1000 of the nearly 80,000 experimental observations. The nature of 
these models for the average omnidirectional counting rate permit 
their closed-form transformation to the equivalent equatorial pitch 
angle distribution. 

The fitted models were sufficiently accurate so that the residuals 
(equal to observed minus fitted) of all the data could be productively 
examined for possible dependence on variables other than the two mag- 
netic coordinates used in the fitting. One consequence of this was the 
detection of instrumental susceptibility to temperature and bias volt- 
age change, which led to an objective partitioning of the data. 

This article summarizes some of the productive aspects of the anal- 
ysis of this body of data. A very large amount of "preliminary" work 
is not reviewed. Though not an historical description of the work, the 
present, paper does have several evolutionary aspects. In particular, a 
series of one-dimensional fits were employed as a base for developing 
two-dimensional models; a preliminary analysis of all the data was 
used to guide the rejection of outliers; a first two-dimensional fit to 
all the data led to a data-independent basis for partitioning the data ; 
the mode of selection of a sample of data, to which the two-dimen- 
sional model was fitted, changed as deeper insight into the importance 
of this issue developed; and, after a very satisfactory fit to the data 
was attained, the model was improved by specialization and reparam- 
eterization so as to overcome some statistical defects and to achieve 
greater physical meaning. 

Readers with specific interests may wish to consult the Table of 
Contents, the summary (Section XIV) and the following overview for 
guidance. 

Section II introduces the input data which have been analyzed. Co- 
ordinates and notation are tabulated, the distribution of the data is 
displayed, and the general quality and stability of the data are dis- 
cussed. It is shown informally that the measurements may be usefully 
organized in the dipole magnetic coordinate system used. 

In Section III, various alternative coordinate systems and scales are 
considered. The bases for choice of the x,L coordinate system for the 
independent variables and the square-root-of-counting-rate scale for 
the dependent variable are given. 

Section IV brings together the ideas underlying the formulation and 
evolution of the models, and gives mathematical definitions and details. 
Some properties of the models which make them suitable smoothing 
functions for this body of data are indicated. 

One-dimensional fits to the data in each of several L-slices (an 
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L-slice is a particular grouping of the data) are displayed on several 
scales and discussed in Section V. It is shown that L-slice fits suffer 
from fundamental deficiencies, in addition to being inconvenient to 
work with. The results of the L-slice fits are used to lead to a two- 
dimensional model. 

Section VI contains the treatment of the preliminary fit of a two- 
dimensional model. This fit is of good quality and provides residuals 
which are used to help identify and eliminate extraneous sources of 
variability in the data and to serve as a basis for more refined sample 
selection. 

The treatment of the two-dimensional fit to the data after it has 
been partitioned to reduce instrumental effects appears in Section VII. 
The method of sample selection is important, and some algorithms and 
their influence on the resultant fits are considered in Section 7.1. The 
advantages of selecting a sample on the basis of a preliminary fit are 
discussed. The fit itself is described and evaluated in the remainder of 
the section. 

A more, detailed statistical critique of the fit discussed in Section 
VII is contained in Section VIII ; in particular, some remaining phys- 
ical and statistical defects are pinpointed. 

Section IX deals with a modified version of the model, which elimi- 
nates the remaining defects, and gives the results of fitting the most 
satisfactorily parameterized model of the proton distribution. 

Residuals are used to study temporal effects in Section X. An in- 
crease in intensity near L = 2 is noted during October, 1962. An upper 
limit of 0.003 gauss is found for the diurnal variation of the earth's 
magnetic field near L = 1.5. A possible shift in the location of the 
atmospheric cutoff is examined. 

The behavior of the radiation belt near the top of the atmosphere is 
the subject of Section XI. Although the data do not allow the position 
of the low-altitude cutoff to be accurately determined, the qualitative 
behavior precludes a simple atmospheric cutoff mechanism. 

Section XII is devoted to a comparison of the Telstar® 1 results, 
presented as flux maps, with those obtained on Injuns 1 and 3, Ex- 
plorers 4 and 15, and other satellites. Absolute flux values agree to 
within a factor of 2 in most cases, which is as well as can be expected. 
Very good agreement exists concerning the behavior of the intensity 
in the equatorial plane, on L-shells, and near the top of the atmos- 
phere. Experimental results regarding the equatorial pitch angle (see 
Fig. 1) distribution are found to agree well with each other, but differ 



PROTON" DATA FROM TELSTAR 1 



1307 



appreciably from the published results of theoretical calculations. 

Section XIII gives brief consideration to possible directions in which 
this work might be extended: improving the fit to the Telstar® 1 high- 
energy protons still further; approaching model development differ- 
ently; employing the data more fully; and encompassing other more 
complex bodies of data. 

Section XIV contains a brief summary of the results and Section XV 
contains acknowledgments. 

Appendix A provides a detailed description of the particle detector 
and its calibration. 

Appendix B gives some statistical background and details of the 
analysis, and Appendix C discusses statistical measures of the good- 
ness of fit of the model over all the partitioned data. 



B = 0.0266 GAUSS 



COS a = 0.5 
a = 60° 




Fig. 1 — Magnetic coordinates of the point P. The spiral is the orbit of a 
particle trapped on the magnetic line of force L = 2.5 and mirroring at B = 
0.0266 gauss. The equatorial pitch angle, «o, is the angle between the velocity 
vector and the magnetic field vector at the equator. 
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II. THE DATA 

The data which are studied in this paper were obtained with a 
detector on the Telstar® 1 satellite which measured protons with energies 
greater than 50 MeV. The sensitive detecting element is a semiconductor 
diode developed specifically for satellite experiments. 7 The effective 
geometric factor, g, of the detector depends upon proton energy, but 
over the region energy between 50 and 130 MeV the average geometric 
factor, g, is relatively insensitive to the energy spectrum and an ap- 
proximate value of 0.143 cm 2 steradian has been selected. These con- 
siderations are discussed in detail in Appendix A. The response of the 
detector is also dependent upon both temperature and electrical bias 
because of changes in the effective thickness of the active region of the 
detector. These effects are discussed in Section 6.8. 

The primary input to our data reduction process consisted of: the 
telemetry record of the number of counts measured by the detector 
in an 11-second counting interval once every minute; the time at which 
the data were recorded (inserted by the recording station); and the 
ephemeris of the satellite position obtained from tracking data. These 
are supplemented by the satellite spin-axis orientation obtained from the 
mirror flash data 8 and by telemetered measurements of the satellite 
skin temperature near the detector and of the detector bias voltage. 

During data reduction, the square root of the counting rate was 
computed for each recorded particle-counting interval and associated 
with the following information: date and time, geographic position, 
position in the earth's magnetic field, orientation of the detector relative 
to the magnetic field, bias voltage, and skin temperature. 

The model developed in the present paper is based on the use of 
a two-dimensional magnetic coordinate system, in which the earth's 
magnetic field is mapped onto an axially symmetric dipole field using 
the adiabatic invariants of particle motion. 9 Any of a number of equiv- 
alent pairs of magnetic coordinates, including the B,L; R,\ and x,L 
sets 10 may be used to locate position in this dipole field. Briefly: The 
magnetic shell parameter, L, specifies a particular line of force (about 
which the trapped particle spirals) by the radial distance to the line 
in the equatorial plane of the dipole measured in units of one earth 
radius (see Fig. 1) ; position along the line of force is specified by either 
the magnetic induction (field strength), B, or by x, where x- (1—B /B)' 
is a convenient variable in the equations of the dynamics of charged 
particle motion. (B is the magnetic induction at the equator on the 
line of force in question.) Magnetic dipole polar coordinates R and X, 
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Fig. 2 — The spatial distribution of data for L < 3 in R, X coordinates. Every 
twentieth point from the //-ordered data is plotted. 



where R is the radial distance in earth radii and X is the latitude angle, 
offer a sufficiently close analog to geographic coordinates to be con- 
venient in many circumstances. The choice among these sets is dis- 
cussed in Section III, as are the reasons for choosing the square root 
of the counting rate as the scale for the dependent variable. 

The coordinates and variables, together with other symbols used 
in this analysis, are listed in Table I under the following headings: 
Radiation Intensity, Position and Orientation, Instrument and Energy 
Spectrum, Mathematical Model, Statistics, and Other. Summary in- 
formation concerning units, constants, derivations, and sources is 
included. 

The satellite was confined to the volume of space {1.09 R, ^ 
R ^ 1.95 R * £ X ^ 5S°j. For \L > 3, R < 1.95 R e \, the average 
counting rate is very nearly zero, and these data were not examined 
further. About 5 percent of the 50-130 MeV proton data for L ^ 3 
were associated with noise bursts which affected adjacent telemetry 
channels; these data were discarded. The study described below is 
based on the remaining 77,649 observations. 

The spatial distribution of the data is indicated in Fig. 2 which is 

*R e = earth radius. 
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a plot in R^. coordinates of the position of every twentieth point from 
the L-ordered data. Although data were not acquired continuously 
during the 226 days that the satellite was active, there arc no time 
gaps in the data longer than two days in duration. 

Fig. 3 is a plot of bands of constant counting rate made by plotting 
the R,\ coordinates at which certain specified numbers of counts were 
recorded during 11-second counting intervals. The data in Fig. 3 cover 
the entire seven-month life of the satellite. The narrowness of the con- 
tour bands demonstrates that the data are exceptionally well-behaved 
in both time and space, and that one may reasonably hope to describe 
radiation intensity in terms of R,k coordinates or their equivalent. 

Among the various sources of error in the data are: noise present 
in the received telemetry signal or introduced during the recording and 
processing of the telemetry; errors in the time as recorded by the 
ground station; errors in the satellite ephcmeris; differences between 
the real magnetic field of the earth and the values of B and L calcu- 
lated from the coefficients in the computer program INVAR (see Table 
I) ; and instrumental effects. In addition, one expects statistical fluc- 
tuations in the measured counting rate at a fixed position. The im- 
portance of these sources of error is discussed later. 




Fie 3 — Bands of constant numbers of counts in 11 seconds in R,\ spare: 
Band a, 4; Band b, 32; Band c, 127-129; Band d, 254-258; Band 3, 508-516 counts. 
All the data from the seven-month period are displayed. 
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in. CHOICE OF THE PRINCIPAL VARIABLES AND THEIR SCALES 

The current state of knowledge of the earth's radiation belts sug- 
gests that the spatial distribution of high-energy protons may reason- 
ably be organized on the basis of a two-dimensional magnetic coordi- 
nate system, except perhaps at very low altitudes near the South 
American magnetic anomaly, where longitude also becomes important. 
Telstai-® 1 data plotted in Fig. 3 indicates that the observed counting- 
rate data does indeed depend principally on the magnetic coordinates, 
R and A. The coordinates R,X are defined in terms of the mathemat- 
ically equivalent pair B,L.° A third equivalent set consists of L to- 
gether with the coordinate x, suggested by Roberts, 10 defined in Table I. 

We have primarily employed the x,L set in this study because of 
the following considerations: In the adiabatic theory, the mirror points 
of particles do not migrate between magnetic shells. 11 Within any shell, 
the coordinate x is approximately linear in A for A < 30°, and thus the 
near-equatorial data is not "crowded" into a small interval of the 
coordinate, as is the case for B. Moreover, we have been able to de- 
velop simple functional representations of the data in terms of x and L. 

The flux of particles is the variable of greatest physical interest for 
comparing the results of different experiments, calculating physical 
effects of the radiation (such as radiation damage to devices in pro- 
posed orbits), deriving an energy spectrum from experimental meas- 
urements, examining the implications of various source and loss mech- 
anisms, etc. However, the flux is not measured directly and requires 
for its calculation knowledge of the energy spectrum of the particles 
and of the energy dependence of the geometric factor of the detector. 
Even in the present circumstances where the conversion is (under the 
assumptions of Appendix A) quite insensitive to these, we prefer to 
carry out the bulk of the data analysis in terms mathematically equiv- 
alent to the directly observed counting rates. 

From among the possible representations of the counting rate in- 
formation (including counting rate, log counting rate, and square root 
of counting rate) the square root of the observed counting rate, Y, has 
been selected as the dependent variable. On the hypothesis that the 
number of counts in a given 11 -second counting interval at any given 
position in space is a random variable with a Poisson distribution, it 
can be shown that the variance of Y is approximately constant, inde- 
pendent of its average value (see Appendix B.2). The least squares 
criterion has been used in all the estimating procedures; that is, coeffi- 
cient estimates have been selected so that the sum of squares of dif- 
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fercnces between observed and fitted values is minimized. The choice 
of the square root scale, Y, as the scale on which to represent the 
counting rate data makes troublesome differential weighting of the 
data in the least squares fitting unnecessary. Similarly, plots of Y 
versus various variables are convenient since the scatter in Y is ap- 
proximately independent of the value of Y. In fact, the square root 
transformation will make the variance of the observation approxi- 
mately independent of its average value whenever the variance is pro- 
portional to the mean. Thus, the procedure is more robust than the 
assumption of a Poisson distribution, for which the variance equals 
the mean. Further discussion and detail is given in Appendices B.2 
and C. 

The results were restored to counting rate and the flux was calculated 
using the best estimate of the average geometric factor, g, (see Appendix 
A) to facilitate the discussion of the physical significance of the meas- 
urements. 

IV. THE EVOLUTION OF THE MODELS 

4.1 General Approach 

This section provides a summary overview of the evolution of the 
models, the details and accomplishments of which are elaborated in 
the following sections and appendices. 

The approach to model development in this study has been largely 
empirical. Theoretical physics considerations are currently too com- 
plex and speculative to do more than serve as a general guide and 
stimulus. We have proceeded on the presumption that an adequate 
model for the spatial distribution of the high-energy protons can be 
based on the mapping of the earth's magnetic field onto a two-dimen- 
sional axially symmetric dipole field, expressed, for example, in the 
coordinates x and L. This is supported by the plots of Fig. 3, the suc- 
cessful polynomial fits on Z,-lines of Mcllwain, 18 Valerio, 10 and Fil- 
lius, 20 and by the results of the present study. 

The ultimate justification of the mathematical models developed 
herein is that, when appropriate estimates of coefficients are inserted, 
good fits to the data are obtained. Various other mathematical, phys- 
ical, and statistical considerations also provided guidance and evalua- 
tion. 

The evolution involved successive interactions with the data and 
iteration on models. Roughly, the main stages included: grouping the 
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data into L-slices; inferring a mathematical function having adjustable 
coefficients which would fit a selected series of L-slices; developing a 
mathematical function to describe the dependence of the L-slice coeffi- 
cients on L; thence fitting the two-dimensional model so-defined to a 
sample of the data; using this fit to screen outliers, to detect instru- 
mental effects and, after partitioning the data, to select a representa- 
tive sample of partitioned data for further fitting; after obtaining a 
very good fit to the partitioned data, some remaining physical and 
statistical defects of the model were overcome by a reparamctrization 
and specialization. Further generalizations of the model were also 
tested. 

4.2 The L-slice Model 

As a developmental operational procedure (encouraged by the L-shell 
orientation of the adiabatic theory 11 ) the data were grouped into a 
series of narrow bands according to L values (e.g., 1.849 ^ L ^ 1.851) 
and plotted versus x. Retrospectively, there is every reason to believe 
that an initial approach based on grouping the data into x-slices would 
also have led to an effective analysis (see Section 13.2). Various func- 
tional forms, having adjustable coefficients dependent on L, were tested 
for adequacy of fit to the L-slices. 

Initially, we employed the functional form 

fc.flfr!«.,fl (*£*,), (1) 



[0 (x > x e ), 

where A, x and S are fitted coefficients for each L-slice, and 

a - *v( i - ifJT' (x = *•>• 

(x > &«)• 

For this body of data from the region {R ^ 1.95 R. , 1.15 ^ L ^ 3.0}, 
we have found this y L {x) function provides an adequately flexible model 
on L-slices, for appropriately fitted values of the coefficients A, x et 
and S. In this representation for given fixed L, the quantity A 2 may be 
interpreted as the average equatorial omnidirectional counting rate, 
since x = on the equator, x e represents a "cutoff" value for x, i.e., 
the cosine of the equatorial pitch angle corresponding to the "loss cone", 
and S has the effect of a shape factor in the y,x dependence. 
The analysis using this y L (x) model is described in Section V. 
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4.3 Dependence on L 

The y L (x) model was fitted to a series of L-slices, obtaining fitted 
values of A, x c and S. These were each plotted against the nominal 
(mid-range) L value for the slice and a reasonably smooth variation 
with L obtained. 

Thence we inferred the following functional dependence of the L- 
slice coefficient estimates on L : 

S = S(L) = s + Si L, (3) 

x e - x e (L) = Vl-ffi'^-siJ 1 , (4) 

R c = R C {L) = L + r,(L - L„) + r 2 (L - L (1 ) 2 + r,(L - L„) 3 , (5) 



A = A'iL) = 






(6) 



10 {L < L ), 

where s , s, , n, r 2 , r 3 , a, , a 2 , a s , •>? and L are fitted coefficients. 

Equation (4) simply expresses the mathematical relationship be- 
tween R (or R c ) and x (or .r ( .) in the magnetic dipole field (sec Table 
I). The coefficient L 0j which occurs in A'(L) and x c (L), may be inter- 
preted as the lower bound of the L shells on which protons with ener- 
gies above 50 MeV were measurable. The quantity R C (L) is such that 
R C (L) — 1 is the equivalent dipole altitude at which the counting rate 
falls to zero. 

4.4 A Two-Dimensional Model — Model I 

The conjunction of (1) to (6) defines a two-dimensional model, re- 
ferred to henceforth as Model I, 

y'{x, L) = A'{L) -G'(x, x e (L), S(L)), (7) 

where G' is essentially the function G of (2) , with x c and S explicitly 
dependent on L. 

Though empirical considerations mainly guided the choice of these 
functions, some physical and mathematical properties influenced the 
choice. In the present case, in which the geometric factor of the de- 
tector is considered to be independent of the energy spectrum (see 
Appendix A), [y(x, L)]- transforms in closed form to the equatorial 
pitch angle distribution, giving 10 
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|2S(L) 



Air 



[A(L)] 2 <1 - 



i 1 - brf} 



fa' L) - g HM 1 + 2S(L)) 



(8) 



where ;(/xo, L) is the predicted equatorial unidirectional flux (protons/ 
cm 2 sec ster) at equatorial pitch angle a = arc cos im, and /3 is the 
beta function, 



/3(p, 5) = f u v ~\l - u)'- 1 du. 

Jo 



(9) 



In addition y'{x, L) has good boundary behavior. The derivative at 
the magnetic equator, d?/(0, L)/dx, is 0, which provides continuity. 
When I < S(L) < J, then dy'(x C) L)/dX -* -oo and d[y'(x c , L)] */ 
dz = 0. The estimated values of S do satisfy this constraint in the 
present case. The desirable consequences of this behavior of the de- 
rivatives will be discussed in Section V. The function y'(x, L) gives 
smooth interpolation over regions sparse in data, and does not have 
any of the wild fluctuations often associated with polynomial fits. 

The analysis of the data using Model I is described in Section VI. 

4.5 Summary Uses of Model I. 

The unspecified coefficients of Model I were estimated by nonlinear 
least squares fitting to a sample of about 1000 observations from the 
complete body of data. Thence this fit of Model I (the CB fit) was 
evaluated relative to all the data and to auxiliary variables, such as 
time, which were not included in the model. Outliers were thereby de- 
tected and screened. An instrumental effect was uncovered (see Section 
6.8), and this led to an objective partitioning of the data, yielding a 
subset (HTB data) for further analysis. The CB fit of Model I was 
also used to specify a representative data sampling procedure for fur- 
ther fitting to the HTB data. 

Though Model I produces a very good fit to the HTB data (see Sec- 
tion VII), it has certain physical and statistical defects. Specifically, 
though the quantities A and x c in the L-slice model have a direct phys- 
ical interpretation, most of the coefficients in y'{x, L) do not. Addi- 
tionally, the estimates of the coefficients in A'{L) turn out to have 
exceedingly high statistical correlations and the model y'{x, L), as a 
function of the coefficients, exhibits marked nonlinearities even in a 
close neighborhood of the least squares estimates (see Section 8.5). 

Therefore, after clarifying the character of the data and obtaining 
a good fit, attention was given to additional improvements of the 
model. 



f 


A,(L - L„) 


l( " - 2) (L - 
,0 


2 [(L p + L- 2L )/2]' 
L,,) + v (L, - L o) - 
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4.6 yi Modified Model— Model II 

The statistical difficulties of Model I were entirely overcome by em- 
ploying a specialized version of A'(L), defined below. Furthermore, 
this specialized model, Model II, retains all the desirable properties 
of Model I while providing both aesthetic improvement and greater 
physical interpretability. 

Model II is defined by 

y"(x, L) = A"(L)-G"(x, x c (L), S(L)), (10) 

where G" is as in (2) , but with S(L) = s , and 

(L ^ L ), 

■!"</■' - " ; " >/,. ^ /-•-! " '-': ' . - : ': -' (id 

(L < L ), 

where A p , L , L p and -q are the coefficients to be estimated. 

A"{L) is a special case of A'(L) and relates to it by the following 
transformations: 

L = L 

v = v 

a 3 = 2L - L P (12) 

a 2 = 2"~\ v - 2)(L P - L )" 

a, = 2'-M„ J? (L„ - L o y~ l . 

Indeed, Model II is essentially defined by the following nonlinear con- 
straint imposed on Model I : 

a 2 = 2"- 1 (, 7 - 2)(L - a 3 )\ (13) 

The coefficients of A"{L) in Model II have the following physical 
interpretations: 

L (as before) is the smallest value of L such that high-energy 

protons are measurable by the instrument; 
A p is the square root of the maximum counting rate of high-energy 

protons in the radiation belt; 
L p is the value of the magnetic shell parameter (on the equator, 

x = 0) at the highest radiation intensity; 
77 may be interpreted as a shape factor for the equatorial (counting 

rate) 1 function, A"(L). 
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The model A"{L) has the form of a product, with the maximum 
value, A p , being multiplied by a factor which decreases as L departs 
from L p in either direction. Note that the factor multiplying A p is 
dimensionless. 

The other fitted coefficients of Model II are s , which is a shape fac- 
tor for the dependence of (counting rate) * on x at constant L, and r,, r 2 
and r 3 which, with L , define the cutoff function x c {L) . 

The analysis of the HTB data using Model II and comparisons of 
Models II and I are considered in Section IX. 

4 .7 Generalization s 

The previously defined models may be regarded as special cases of 
Model III defined by 

y'"(x, L) = A'"(L)-G'"(x, x£L), M(L), P(L), Q(L)), (14) 
where A"\L) = A'(L), defined in (6) , 



b - (£)' 



P(L) 

2\ 0(1.) 



/(i - x>r L) (x ^ x.), 



(15) 



.1/ ( /. ) 

! i f -^-1 

G'" = - 

[0 (* > x c ), 

x r (L) is as defined in (4), and M(P), P(L) and Q{L) involve coef- 
ficients or functions to be fitted. 

The function G' is a special case of G" f , in which M(L) = 2 and 
Q(L) = $. This permits a closed form transformation to an equatorial 
pitch angle distribution. The function G" additionally constrains P(L) 
= s , independent of L. 

The more general G'" in Model III can be used on L slices to de- 
termine L-slice estimates of M, P, Q, as well as ^1 and x c , and these 
in turn inspected to infer functional dependence on L. Clearly, this 
more general form must lead to at least as good a fit as Models I or 
II. Work has been done with Model III 21 but no important improve- 
ment over Model II was obtained for this body of data. 

Neither of the fitted models y'(x, L) nor y"{x, L) is applicable far 
outside the spatial and energy regions that include the data analyzed 
here. For example, Models I and II do not fit well to the 26-33 MeV 
protons measured by the Telstai-® 1 satellite, nor are they suitable for 
fitting many of the electron distributions. Preliminary investigations 
indicate that these remarks may not apply to G'", whose additional 
coefficients allow more rapid changes in curvature as a function of x. 
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We have already shown for Telstat® 2 data 5 that A(L) can be ex- 
tended to include description of the plateau of high-energy protons 
reported by Mcllwain 18 - " near the equator at R S 2.2 R , beyond the 
orbital extremes of the Telstar® 1 satellite. The extension was made by 
adding a term to A'(L), (6), to give A lv defined by 



A" = A'(L) + a 4 exp 



sk^q , (16) 

a 5 J 



where a 4 , as, and L\ are coefficients describing the equatorial distribu- 
tion of the "excess" protons that give rise to the plateau. In the less 
stable parts of the radiation belts the early work on empirical time 
dependence presented by Gabbe and Brown 5 clearly requires extension. 

V. FITS ON THE L-SLICES 

The model of (1) and (2) was fitted to the data, on the scale of Y, 
in 92 individual L-slices, using a nonlinear, multidimensional, least 
squares, computer program (see Appendix B) to estimate the coeffi- 
cients and produce various statistical measures. The procedure of fit- 
ting to L-slice data enabled one to test functional forms of j/j [x) and 
then to evolve functional forms for the dependencies of the coefficients 
of the L-slice models on L. 

Proceeding in this manner, however, has a number of possible pit- 
falls. In particular, the estimates of coefficients within an L-slice may 
be highly correlated, and the reliability of the actual values of the 
estimated coefficients also depends on the pattern of data points in 
the particular //-slice, e.g., whether or not there are points near x c . 
Hence, the estimated values for any particular coefficient may not ex- 
hibit a smooth dependence on L. 

The form of the L-slices whose middle values of L, called L,„, are 
1.35, 1.801, 2.2015, and 1.79, respectively, are displayed in Figs. 4 to 7. 
The thin solid lines in the figures are the fits to the 7,-slicc data (mean- 
ing of the dashed and thick solid lines will be taken up later). The 
numerical values of the coefficients of the fits, and the widths of the 
slices are given in Table II. Figs. 4 and 5 are examples of the high 
quality of fit which is typically obtained for L-slices having L,„ < 2. 

In Figs. 4(a) and 5(a), square root of counting rate is plotted 
against x. One sees that the fit to the data points (the thin solid line) 
is quite adequate. The cutoffs, .r,., are well-defined, the scatter in Y is 
approximately independent of y and the data arc well-distributed in x. 
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Fig. 4 — Data from the L-slice centered at L m = 1.35 and the results of three 
fits shown on four scales. 
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Fig. 5 -Data from the L-slice centered at L m = 1801 and the results of three 
fits shown on four scales. The partitioning in (a) is discussed in bection 7.1. 
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Table II— Coefficients and Statistics of the L-Slice Fits. 



Lm 


1.35 


1.801 


2.2015 


1.79 




1.346 


1.800 


2.200 


1.7895 






1.354 


1.802 


2.203 


1.7905 




AL 


0.008 


0.002 


0.003 


0.001 




A 


6.757 


4.109 


1.70 


4.324 




<r(A) 


0.053 


0.031 


0.12 


0.043 






0.6795 


0.8998 


0.954 


0.923 




<r(x c ) 


0.0027 


0.0044 


0.011 


0.015 




S 


0.324 


0.390 


0.58 


0.478 




t(S) 


0.018 


0.024 


0.10 


0.060 




Number of pts 


140 


129 


144 


65 




MSE 


0.1125 


0.0497 


0.0282 


0.0478 




Correlation coefficients 


A with x c 


0.281 


0.309 


0.724 


0.408 




A with 5 


0.605 


0.561 


0.940 


0.548 




x c with S 


0.774 


0.820 


0.890 


0.944 





As the cutoff is sharp on the scale of y, it is convenient to have a 
function which has an infinite derivative at x„. Otherwise the exact 
x at which y -» may have relatively little effect on the mean square 
error of the fit. This would lead to an ill-defined value for x c , even 
though the data allows one to evaluate the position of the cutoff quite 
precisely for L values smaller than si. 9. 

In Figs. 4(b) and 5(b), the counting rate, Y 2 , is plotted against x. 
The thin solid lines represent the same fits as those in Figs. 4(a) and 
5(a). One finds that the position of the cutoff is no longer well-defined 
on the plot. Instead the counting rate fades away as x increases. Hav- 
ing the derivative of y 2 equal zero at the cutoff (as noted in the pre- 
vious section) is suitable in this situation. The scatter in Y 2 now 
changes with y 2 , and is greater for large values of y- (small values 
of x). This nonuniform scatter makes it more difficult to judge the ap- 
propriateness of fit. If one wished to minimize the squared deviations 
between observed and fitted in terms of y 2 (or log y 2 ) the values of 
Y 2 (or log Y 2 ) would have to be weighted inversely as their estimated 
approximate variance, with a loss of intuitive appreciation of the qual- 
ity of fit from a scatter plot and a substantial inconvenience in carry- 
ing out the fitting procedure. 

In Figs. 4(c) and 5(c) the ordinate is log y 2 . This choice of coordi- 
nate restores the ability to discriminate in the vicinity of the cutoff at 
the cost of a large loss of sensitivity in regions where the counting rate 
is higher. 

Finally, Figs. 4(d) and 5(d) display the same data in the coordinate 
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system log y-, log (B/B ) . This choice of abscissa expands the high-.r 
region enormously, but contracts the low-.t region to the point where 
it is impossible to see the details of the particle distribution in the 
vicinity of the equator (x = 0). This contraction would be even more 
severe if the abscissa were B or B/B„. 

In the region defined by A < 45°, which covers the high energy pro- 
ton data, the coordinate x provides adequate detail (see Ref. 10 for 
further discussion). If, however, the data had extended to A > 45° an- 
other choice of magnetic coordinate would have been desirable for 
X > 0.95, because all A > 45° are crowded into x values between 0.95 
and 1. 

The standard errors and correlations of the coefficients of the 
four L-slices under discussion, together with mean square error (MSE) * 
of fits, are listed in Table II. The standard error is in general a 
relatively small fraction of the estimate and the MSE is substantially 
greater at small values of L m than at larger ones. This is further ana- 
lyzed in Section VI. 

At L = 2.2 the satellite gets no closer to the magnetic dipole equator 
than A = 20°. This fact, which is associated with the problem of cor- 
relation of coefficient estimates within L-slices, is displayed more em- 
phatically by choosing x as a coordinate, as in Figs. 6(a), (b), and (c), 
than by choosing log (B/B Q ) as in Fig. 6(d). In addition, in Fig. 6(d) 
the expansion of the abscissa in the region of the cutoff makes it diffi- 
cult to judge the physical appropriateness of the value of x which re- 
sults from the least squares procedure. The same difficulty is encoun- 
tered to a lesser degree with Fig. 6(b). However, in Figs. 6(a) and 
6(c) one judges the .r-intercept of the thin solid line to be too large, 
and Fig. 6(a) has the additional advantage of allowing one to make 
a better judgment of the quality of the fit at lower values of x. As might 
be surmised from the high values of the correlations for L m — 2.2 in 
Table II, the value of x can be adjusted to a substantial extent with- 
out much change in the mean square error. These high correlations, 
which typically occur for L, n > 2, reduce confidence in the individual 
estimates of the coefficients for given L-slices. This difficulty also re- 
duces the stability of the estimates of the coefficients as L m is changed, 
and precludes basing the values of x c {L) and S(L), for L > 2, on the 
fits to the L-slices. 

A similar difficulty may be introduced when L < 2 by sampling 
fluctuations as illustrated in Fig. 7. In this case, there is a scarcity of 



* Some statistical terms are defined in Table I. 
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Fig. 6 — Data from the L-slice centered at L« = 2.202 and the results of three 
fits shown on four scales. 
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Fig. 7 — Data from the L-slice centered at L m = 1.790 and the results of throe 
fits shown on four scales. 
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data near and beyond the cutoff, unlike the slice with L m = 1.801 
illustrated in Fig. 5. The paucity of data near the cutoff in the L-slice 
centered on L m = 1.79 both correlates and distorts the values of .r, 
and S. In this particular case, the width of the L-slice can be increased 
to avoid this difficulty, but, in general, increasing the width of the 
slice to include enough data may introduce a serious L-dependence 
within the slice. As a result, x may be determined by points near one 
extreme of L within the slice, .4 by points at the other extreme and S 
by some combination. This problem is especially severe below L = 1.3 
where data begin to become sparse. 

The plotted points in Figs. 8 to 10 summarize the dependencies of 
the estimates of the L-slice coefficients A, x , and S, respectively, on 
L M , for all 92 slices. More than one value of the coefficients is plotted 
for some values of L,„ because on occasion the width of the L-slice was 
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Fig. 8 — Three estimates of A as a function of L. 



PROTON' DATA FROM TELSTAR 1 



133" 



08 



06 




O L- SLICE ESTIMATES 

FROM CB COEFFICIENTS 

FROM HTB COEFFICIENTS 



Fig. 9 — Three estimates of x c as a function of L. 



varied without changing L m . Although there are local fluctuations in 
the estimates that arise from the way a narrow L-slice samples the 
data, the estimates exhibit a smooth dependence on L. The fluctuations 
are particularly pronounced near L m = 1.8 in Figs. 9 and 10, and L„, = 
1.3 in Fig. 10. 

The standard errors of the L-slice estimates of A are typically 1 per- 
cent for L < 1.95, but become as large as 6 percent where there are 
no equatorial data, as is the case for L > 1.95. Fox x e estimates, the 
standard errors are typically 0.5 percent. The estimates of S have a 
standard error of about 5 percent (±0.015) near L = 1.5 and about 
15 percent (±0.05) near L = 1.2 and L = 2. The meanings of the 
curves in Figs. S to 10 will be discussed in the following sections. 
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Fig. 10 — Three estimates of S as a function of L. 

In summary, the. L-slice approach enables one to infer a functional 
dependence of L-slice coefficients on L and to obtain an intuitive ap- 
preciation of the quality and nature of fit. The fitting procedure re- 
quires refinement by being carried out as a simultaneous two-dimen- 
sional process in x and L jointly. This overcomes the "grouping" 
inaccuracy in the L-slice approach and in addition makes good use of 
the data in those regions where data are scarce. The resultant function 
also provides convenient and excellent interpolation of data over the 
entire x,L region while employing a relatively small number (8, 9, or 
10) of fitted coefficients. 



VI. THE TWO-DIMENSIONAL FIT FOR THE COMPLETE BODY OF DATA 

The analysis of this section is a precursor to the more refined paral- 
lel analysis of Section VII. This preliminary analysis produces the 
following results of consequence: Model I (see Section 4.4) is shown 
to be satisfactory; instrumental effects are identified and an objective 
algorithm for partitioning the data to reduce these effects is formu- 
lated; outliers are screened; and a more adequate basis for sample 
selection is provided. Many statistical details are omitted from this 
section, and statistical matters are dealt with more fully in Sections 

VII, VIII, and IX and in Appendices B and C. 

6.1 Sample Selection and Fit 

It was necessary, for practical computing reasons, to make a selec- 
tion of approximately 1000 observations on which to carry out the 
simultaneous two-dimensional (in x and L) nonlinear (in the coeffi- 
cients) least squares fit. In this preliminary phase, the nearly 80,000 
data points were sampled by dividing the L-range from 1.15 to 3.00 
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into 925 contiguous intervals, each 0.002 wide. One data point was 
selected from each interval. As the data are approximately uniformly 
distributed in .r (in the .r-range covered by the satellite) in each 
L-slice (see. Figs. 4 to 7), no effort was made at this point to in- 
fluence the x distribution of the observations in this subset. The ques- 
tion of the "design" of the sample to be used as a basis for fitting the 
model is rather important, however, since the fit obtained with the 
empirical model is responsive to the distribution of data in x,L space. 
Other bases of sampling were employed later (see Section 7.1 and Ap- 
pendix B.3). 

Model I, described in Section 4.4, was fitted to the 925-point sample 
from the complete body (CB) of data. As this serves only as a pre- 
liminary fit, the values of the CB coefficients and other statistics are 
not presented here. 

The quality of this fit was examined from various viewpoints: (i) 
by its behavior along the boundaries of the belt; (ii) by comparison 
with the //-slice fits; (Hi) by plotting the residuals (observed value 
minus fitted value) versus the x and L coordinates; and (iv) by ex- 
amining the mean square residuals (MSR) in various regions of mag- 
netic coordinate space. Though the coefficients of the model were esti- 
mated from 925 sampled data points, the evaluation of quality of fit 
was based on all the nearly 80,000 observations. 

6.2 Evaluation of Fit at Equator 

The points in Fig. 11 are the values of Y (square root of observed 
counting rate) plotted against L for all data points for which x is near 
0, specifically x < 0.037 (i.e., A. < 1°). For a given L, y'(x, L) changes 
very little between x = and x = 0.037 (see Figs. 4 and 5) and the 
points in Fig. 11 may be regarded as approximate equatorial points. 
The curve in Fig. 11 gives the fitted values of A'(L) = y' (0, L) using 
the CB coefficients, and appears to represent the data very well. Note 
that A'{L) has not come from a fit to the equatorial data as such, but 
rather is the equatorial value of if as predicted by the two-dimensional 
fit. That is, the fitted A'{L) does not minimize the sums of squares of 
deviations for just the equatorial points, but is, rather, the optimum 
fit in the least squares sense to the 925-observation sample, and these 
observations are distributed through x,L space. The excessive scatter 
in the equatorial value of Y between L = 1.35 and L = 1.55 which 
shows in Fig. 11 will be taken up in the next section. 

The values of A'iL) are also plotted for reference as the dashed 
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Fig. 11 — All data for x < 0.037 (i.e., within 1° of the magnetic invariant 
equator) and the equatorial value estimated from the CB coefficients plotted 
against L. A' and Y are in units of (counts/sec) 1/a . 

line in Fig. 8. One sees that the L-slices give quite good estimates for 
A, although these estimates tend to be a little erratic and to favor 
the lower values rather too much in the neighborhood of L = 1.4. 



6.3 Evaluation o] Fit at Cutoff 

The cutoff may be thought of as the position of the outer envelope of 
the nonzero counting rate, or the inner envelope of the zero counting 
rate. Thus, in practice the location of the cutoff is associated with the 
sensitivity of the detector, rather than with the absence of particles. 
For L ^ 2, there is a wide range of x over which there are many in- 
stances of either zero or one count occurring during the 11-second count- 
ing interval, and as a result the cutoff is not well-defined. This is 
exemplified in Fig. 6. The overlapping of the region in which no count is 
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observed with that in which one count is observed shows clearly in 
Fig. 12. The locations of occurrences of zero counts are plotted in R f \ 
coordinates in Fig. 12(b) and in x,L coordinates in Fig. 12(d). Figs. 
12(a) and (c) show the locations at which one count (one, two, and three 
counts for L < 1.5) was recorded. (The density of points has been re- 
duced at high L to improve the clarity of the display.) 

Because the cutoff is increasingly difficult to define from the data as 
L increases beyond p^2, the position of the cutoff predicted by the fitted 
model is not a good boundary condition to use in judging the quality of 
the two-dimensional fit. Instead the locus of positions for which exactly 
one count per counting interval is predicted is superimposed as the solid 
lines in Figs. 12(a) and (c) upon the array of points giving the band 
of positions at which one count per counting interval was observed. The 
data are represented quite satisfactorily by the solid lines particularly in 
the region (L ^ 1.90) where the belt ends abruptly. The fit is least 
satisfactory near L = 2 (X = 40°). Adding the terms r 4 (L — L ) 4 and 
r 6 (L — L ) 5 to the expansion for R e (L) in (5) does not appreciably 
improve the fit near A = 40°. 

The line x c (L), representing the cutoff itself, is plotted as the dashed 
line in Fig. 12 and is seen to be a reasonable outer envelope for the 
nonzero counts. 

The present estimate of .r,.(L) is also shown as the dashed line in 
Fig. 9. Below L JS 1.8, the estimates of x c from the individual L-slices 
are in good agreement with estimates from the two-dimensional fit. 
However, above L S 1.8 the L-slices give erratic values for x . As 
demonstrated in Fig. 7, the L-slice estimates may be biased toward 
high values, a circumstance which makes it difficult to extract a satis- 
factory fit for x c {L) from the estimates of x c produced by fitting the 
L-slices. 

6.4 Behavior of S(L) 

The values of the function S(L) generated by the two-dimensional 
fit cannot be subjected to a simple boundaiy comparison with the data. 
The function S{L) is plotted as the dashed line in Fig. 10 along with 
the L-slice estimates. It will be seen that the L-slice estimates tend to be 
somewhat higher than the values given by S{L) in the neighborhoods 
of L = 1.3 and L = 1.9. However, if the form of S{L) is taken to 
provide a better fit to the points in Fig. 10, then the resulting two- 
dimensional fit yields a physically less satisfactory fit of the cutoff 
function x c {L) to the boundary data without substantial improve- 
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Fig. 12 — All positions in R,\ space (a) and x,L space (c) at which one count 
(one, two, and three counts for L < 1.5) was observed in an 11-second counting 
interval, and all positions in R, \ (b) and x, L space (d) at which zero counts 
were observed in an 11-second counting interval. The solid lines are the loci of 
positions at which the CB coefficients estimate one count in 11 seconds. The 
dashed lines are the loci of the cutoff function x e (L) or R C (L) calculated from 
the CB coefficients. The trace R = 2.0 Re, which explains the absence of data 
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Fig. 12 — (continued) 

in the lower right-hand corner of the x,L plots, appears in part (d). The cluster 
of points near R = 1.1 and \ = 20° in part (b) of the figure is data acquired by 
the telemetry station at Woomera, Australia. It represents observations made 
near perigee when the satellite was below the bottom edge of the proton belt 
which is high over the western Pacific Ocean. 
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ment in the overall fit (see also Section 4.7). Admittedly, this judg- 
ment is subjective because it is made in regard to regions where the 
cutoff is poorly defined by the data because of the insufficent sensi- 
tivity of the detector. The high values of S near L = 1.9 appear to 
arise from the correlation problem discussed in Section V in connection 
with Fig. 6 and Table II. 

6.5 Behavior of the Fit on Several L Slices 

The dashed lines in Figs. 4 to 7 are the values predicted by the CB 
coefficients superimposed on the L-slice data along with the pre- 
viously derived L-slice fit. In Figs. 4 and 5, the difference between the 
thin solid and the dashed lines is insignificant, and this is generally 
the case for L < 1.95. At L„, = 1.79, the predictions from the CB 
coefficients differ importantly from the fit to the L-slice only for x 
values at which there are no data. 

For L,„ = 2.2, however, the two predictions are noticeably different 
as may be seen in Fig. 6. The fit to the L-slice gives the estimate 
x c = 0.954 (see Table II) ; the two-dimensional fit yields x B = 0.928; 
and the difference exceeds two standard deviations. The question as to 
which of the two lines is a better representation of the data in this 
L-slice in the physical sense, rather than in the least squares sense 
applied to these points by themselves, is connected with criteria 
which will be discussed in the following sections. The basic fact is 
that the two-dimensional fit provides a mechanism by which the data 
on every L-slice can influence the fit on every other L-slice and 
thereby provides a fit that is more satisfactory overall than the 
collection of individual L-slice fits. 

6.G Residuals in x,L Space 

The data were also examined for dependencies on .r. and L over 
and above those provided for by the fitted mathematical model. This 
is accomplished by studying the residuals, i.e., (Y — y), for all the 
nearly S0.000 observations. The residuals provide a very sensitive basis 
for judging the quality of the fit. The removal of the principal depen- 
dence on x and L by subtracting the fitted function from the observa- 
tions has the effect of allowing small systematic differences to be 
prominently displayed. 

Fig. 13 shows a 3100-point sample of the residuals, Y — y, plotted 
against L, where, to keep the density of the points reasonable, only 
one point has been plotted from each of the nearly 3100 contiguous 
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Fig. 13 — CB residuals of F(i.c, Y — y calculated from ,t he CB coefficients) plotted 
against L. The arrows indicate ± the approximate standard deviation if Y- were 
Poisson distributed. No more than one point, is plotted for an L increment of 0.0006. 

L-intervals, of width AL = 0.0006, between L = 1.15 and L = 3. 
Ideally, the residuals should scatter randomly about 0, without any 
perceivable pattern. For L < 2.4 there is only a little indication of a 
no n random trend. However, for L > 2.4 there is a distinct pattern. 
This pattern is associated with the quantization error, which becomes 
important where the number of counts per counting interval is very 
small. When < y < Vl count/ 11 sec and Y = or Vl count/11 sec, 
the result is the tailing upward toward the resi dual = axis t hat starts 
at L tt 2.4. When y = and Y = or yl count/11 sec, one gets 
the two-line pattern (0 and 0.0310 = Vl/ll) seen clearly in Fig. 13 
for L ,> 2.7. (The thickening of the zero axis indicates the presence 
of data points.) 

Fig. 14 is a plot of the residuals against x for all points for which 
1.4 < L < 1.6. The residuals in Fig. 14 show no structure; however, 
their average value is a little less than zero. This dip is confirmed by 
the points in the range 1.4 < L < 1.6 in Fig. 13, and means that the 
value of y is slightly high relative to the data in this region. However, 
the lack of structure in Fig. 14 indicates that the bias is independent 
of x in this region. 

Fig. 15, the plot of the residuals vs x for 1.85 < L < 1.90, shows 
the region in which the fit is poorest. The residual points are not sym- 
metrically distributed about zero and the asymmetry seems to depend 
on x. Notice that the value of y is slightly too large near z ~ 0.05 and 
x ~ 0.65. The discussion of these trends is continued below, after 
some further analysis has been described. 

6.7 Mea?i Square Residuals in x,L Space 

Another way of gauging the finality of fit is to compute the mean 
square of the residuals (MSR) separately for various regions of 
x,L space. Trends in these quantities may indicate regional varia- 
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Fig. 14 — CB residuals of 7 (i.e., Y—y calculated from the CB coefficients) plotted 
against x for 1.40 < L < 1.60. The arrows indicate ± the approximate standard 
deviation if Y 2 were Poisson distributed. 

tions in the adequacy of fit. The data and residuals were divided into 
three groups. Group I contains all the "good" data points "within" 
the boundaries of the > 50 MeV proton belt. These points are defined 
as those not included in Groups II and III. Group II consists of the 
"good" data points "outside" the boundaries of the belt. These are 
points which meet two criteria: they have values of (x, L) for which 
x is greater than x c (L) + 0.001, and they are not in Group III. Group 
III comprises the outliers or "bad" data points, defined as those points 
whose residuals are greater than three times the overall root mean 
square residual of the points in all three groups together.* The most 
probable origin of a point in Group III is a telemetry error. 

If the number of counts in a counting interval behaves like a 



* Note that only 0.5 percent of the data fall in Group III. 
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Poisson random variable, then the variance of Y- would be equal to 
the average value of Y 2 . As noted in Appendix B, when Y is not near 
zero, the variance of Y would then approximately equal 0.023, inde- 
pendent of the average value of Y. This value then might approx- 
imately represent the average value of the mean square residual, 
MSR, on the scale of Y. Thus, the number 0.023 provides a baseline 
for the comparisons discussed below. 

Table III lists the mean square residuals (MSR) by L range and 
by Group. For Group II, }' is frequently zero and, as x > x c implies 
y = 0, one finds that the residual is zero very often. Of course, under 
the Poisson assumption the variance of Y when its average value is 
or very close to will be less than 0.023 (see Appendix B.2) and the 
appearance of MSR values smaller than 0.023 in Group II is thus not 
surprising. A similar circumstance exists in Group I for L > 2.6. 
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Fig. 15 — CB residuals of Y (i.e., Y— y calculated from the CB coefficients) plotted 
against x for 1.85 < L < 1.5)0. The arrows indicate ± the approximate standard 
deviation if I'- were Poisson distributed. 
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For the overall fit, the MSR of Group I (L range from 1.1 to 3.0) is 
only twice 0.023. However, for 1.3 < L < 1.6 the Group-I MSR is four 
times 0.023. This L range is associated with the large scatter in the 
equatorial data plotted in Fig. 11, and Fig. 14 shows that this scatter is 
independent of x, rather than just an equatorial phenomenon. This 
issue is pursued further below. 

6.8 Dependence of Residuals on Other Variables 

Studies were made of the possible dependence of the residuals on 
observed variables other than x and L. Indeed, it will appear that 
some of the excess scatter exhibited in Table III and in Figs. 11 and 
14 is associated with instrumental effects. 

The regularities inherent in the orbit and orientation of a satel- 
lite, the motion of the earth, and the location and operation of the 
telemetry receiving stations lead to systematic interrelations among 
the various coordinates listed in Table I. A simple example concerns 
temperature. The satellite cools when its enters the earth's shadow. 
This eclipse occurs only on the night side of the earth. Thus, if the 
detector is temperature sensitive, one would see a false day-night 
effect in the counting rate. If, because of additional dependencies, 
data are available during eclipse for only a limited span of days, a 
false secular effect might also be observed. Because of the implications 
of the preceding discussion, a careful study was made of the behavior 
of the residuals with respect to a large number of coordinates, and 
attention was given to the details of the relationships among the 
coordinates during the search for contributors to the inflation of the 
MSR. 

We present below the evidence that has led us to the conclusion that 
two instrumental effects, variations in bias voltage and changes in 
temperature of the detector, are principal causes of inflation of the 
MSR. 

There was no temperature sensor on the particle detector. The 
instrument is not exposed to sunlight and is relatively well-insulated 
thermally from the skin and frame of the satellite. Consequently, 
temperature measurements of the skin are not closely related to the 
temperature of the detector. However, a good indicator of detector 
temperature is elapsed time since entering or since leaving eclipse. 
Fig. 16 gives plots of the residuals, Y — y, against time in minutes 
measured from the more recent of the two events, entered shadow or 
entered sunlight. Residuals associated with periods during which the 
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satellite did not enter eclipse once per orbit are segregated at the far 
right-hand side of the plots, labeled .4 on the abscissa. 

Figs. 16(a) and (b) are for 1.4 < L < 1.6. The points in Fig. 16(a) 
are those for which the bias voltage was between 95.3 and 97.5 volts, 
while Fig. 16(b) contains those for bias voltages between 92.0 and 
95.3 volts. The decrease in the residuals (and also in the observed count- 
ing rate) after the satellite enters eclipse (and the temperature falls) 
and the increase after the satellite leaves eclipse (and the temperature 
rises) may be seen distinctly in both figures. In addition the residuals are 
noticeably more negative for the low (92.0 to 95.3 V) bias range. Both 
low bias voltage and low temperature are known to decrease the ef- 
ficiency of the detector and one expects an appreciable effect to be intro- 
duced into the counting-rate data. In the present case the scatter is 
about ±15 percent of the counting rate. A consequence of this is the 
excess scatter that has been noted particularly with reference to Fig. 11 
and Table III. 

Figs. 16(c) and (d) are analogous to Figs. 16(a) and (b), but the 
residuals are for the L range 1.85 to 1.90. Again, the systematic 
influence of low temperatures and low bias voltages is unmistakable. 

6.9 Partitioning the Data 

Two ways of responding to these instrumental effects might be: 
(i) to try to correct the data, or (ii) to disregard the affected data. 
It is not possible to make a correction to the counting rate that is 
properly independent of the experimental results because; (i) the bias 
voltage was measured in steps of 1.11 Y, which is not sufficiently fine- 
grained; (ii) it would be necessary to estimate the temperature of 
the instrument using a complicated hypothetical relationship between 
the instrumental temperature, skin temperature, and time after enter- 
ing eclipse (or sunlight) ; and (Hi) we have an insufficient knowledge 
of the temperature and bias-voltage sensitivity of the detector. 

Though an ad hoc correction based on the observed counting rates 
could have been attempted, it was decided for practical reasons to 
eliminate both the low-temperature and low-bias points and use only 
that data which was gathered under the following conditions: 

(i) The satellite had been in sunlight for the previous 50 minutes, 
and thus had attained temperature equilibrium reasonably well 
(see Fig. 16). 

(ii) The bias voltage was between 95.3 and 97.5 volts. 
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Fig. 16 — CB residuals of Y (i.e., Y — y calculated from the CB coefficients) plotted 
against time in minutes from the most recent of the two events, entered eclipse 
or entered sunlight. Data taken on days during which no eclipse occurred are plotted 
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Fig. 16 — (continued) 

within the region marked "A" at arbitrary values of the abscissa. The arrows 
indicate ± the approximate standard deviation if Y* were Poisson distributed. 
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This selection yields a homogeneous body of 41,135 points, hence- 
forth referred to as high temperature-high bias (HTB) observations. 
The remaining 36,500 points, which represent a mixture of tempera- 
ture and bias conditions, were used only occasionally in further 
analyses. This selection process coincidentally produces one unfort- 
unate associated circumstance, namely, the exclusion, as low-bias 
data, of all measurements made between days 325 and 373. 

Further analysis and model fitting and development based on, and 
directed towards, this HTB data is detailed in the following sections 
and Appendix C. 

VII. THE TWO-DIMENSIONAL FIT FOR THE SELECTED (HTB) DATA 

7.1 Sample Selection 

The distribution of the HTB data in magnetic space is indicated 
in Fig. 17, which gives the #,A coordinates of every tenth point from 
the 41,135 L-ordered HTB observations. The data provide reasonably 
adequate, though uneven, coverage. As a practical requirement for the 
fitting procedure, a "representative" sample of about 1000 observa- 
tions must be selected. 




Fig. 17 — The spatial distribution of the HTB data for L < 3 in B, X coordi- 
nates. Every tenth point from the L-ordered data is plotted. 
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It is intuitively clear from preliminary knowledge of the radiation 
distribution that some sample configurations will be far more effective 
than others in defining the functional form of the proton flux. 

The sample selection is important because: (i) nothing more than 
a sophisticated smoothing function is being fitted and we want this 
function to be broadly applicable over the entire space; (ii) an 
optimum fit in one region of space does not necessarily imply a good 
fit elsewhere; (Hi) the spatial distribution of data points depends on 
the satellite orbit and the position of the telemetry stations; (iv) even 
with the square root transformation, there remains some differential 
variance among the data. 

These considerations argue against using a simple random sample 
or even a random sample in x with a systematic sample in L such 
as in the CB fit. Indeed, they also argue against fitting all (un- 
weighted) HTB data, even if this were practical. Alternatively, points 
might be chosen on the basis of a simple geometric grid in magnetic 
space. Such a procedure would be easy to use, but it is arbitrary with 
respect to the radiation belts. 

Sampling procedures might be based on particular physical features 
of the radiation belts to emphasize the goodness of fit, for example, 
where the flux is high or where diffusion across L lines might be 
important. However, such fits would be too biased for our present 
general objective. 

One is thus led to a sampling process based on properties of the 
radiation belt itself, as described for example by the preliminary CB 
fit. In particular, a high density of data points is desirable in regions 
where the value of y is changing rapidly, while a low density will 
suffice where the function is changing slowly. A realization of this 
criterion would be to define about 1000 x,L cells, within each of 
which the range of y from the preliminary fit would be the same. 
However, there arc appreciable practical difficulties in defining the 
boundaries of such cells. 

Thus, the following hybrid procedure was used to define the 960- 
point HTB sample on which the subsequent fitting was done: The 
L-rangc from 1 to 3 was divided into about 120 L-slices of equal 
(~ 0.017) width in L. Each L-slice was then divided into eleven 
x,L cells using a scheme that depends on the preliminary fit. The 
first ten cells were chosen so that within each cell the range of y 
predicted by the CB model is closely 1/10 of the equatorial value of 
II at the center of the L-slice. The eleventh cell lies beyond x c . The 
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method of partitioning in the x direction is illustrated by the partition 
of the L-slice in Fig. 5(a) into five rc-regions by the horizontal 
lines. (The distance d is added to x to define the lower-.r boundary of 
the last cell.) 

To take some account of differential variances remaining after 
the square root transformation, the following procedure was em- 
ployed: The mean square deviation from the mean (MSD*) was 
calculated for all the HTB data in each x,L cell defined above; 
thence, after visual inspection of the results (see Appendix C) , three 
groupings of contiguous x,L cells were made according to whether the 
MSD's were generally below 0.013, between 0.013 and 0.020, or above 
0.020; the corresponding regions were then given relative weights of 
2, 1£, and 1, respectively. The weight 1 implies that one point was 
sampled from the cell. 

These weights were assigned on the basis of a judgment which con- 
sidered: (i) the desire to increase the weight of low variance (i.e., 
near-zero counting rate) observations and thus to aid the definition of 
the cutoff; and (ii) the desire to keep from "wasting" sample points 
in the region x » x since such data will add little to the specification 
of x c (L) and virtually nothing to the estimation of A{L) and S. 

Fig. 18 shows the distribution in x,L space of the 960-point sample 
which was used. The number 960 came about because a number of 
the defined cells had no data in them. Our experience with several 
other samples of the HTB data gives us confidence in both the ration- 
ale behind, and the results obtained with, this 960-point set, henceforth 
referred to as the HTB sample. However, sampling procedures tailored 
to the requirements of special purpose fits will give better results in 
some regions of x,L space. 

Some additional discussions relevant to sample selection and data 
usage are given in Section 13.3 and Appendices B.3 and C.2. 

7.2 The HTB Fit 

A slightly constrained version of Model I of Section 4.4 was fitted 
to the 960-point HTB sample. The results are referred to as the HTB 
fit. The constraint is Sx = 0, in (3). Most of the values of Sj. obtained 
in preliminary fits to various samples of the HTB data differed from 
zero by less than two standard deviations. Also, the points in Fig. 



* See Table I for definition of MSD, MSR and MSE. 



PROTON DATA FROM TELSTAlt 1 



1357 




Fig. 18 — The distribution of the 960-point HTB sample in x, L space. The 
trace It = 2.0 fl« explains the absence of data in the lower right-hand corner of 
the figure. 

10 do not suggest a linear dependence of S on L* The effect of this 
constraint on the value of the fitted cutoff function was examined and 
found to be unimportant. 

The estimated HTB coefficients (obtained by fitting the constrained 
model to the HTB sample) appear in Table IV. The physical inter- 
pretation of L as the lowest L on which > 50 MeV protons were 
measurable was noted in Section 4.3. The standard error of 0.001 (~6 
km in altitude) is no larger than the uncertainties inherent in the 
calculation of L itself. 

The interpretation of S as a shape factor (see Section 4.2) is 
straightforward in the present case, i.e., where Sj = 0. The standard 
error of 0.005 is much smaller than the standard errors of the 
estimates of S generated from the fits to L-slices (Table IT) and is 



* Some higher-order models for S(L) were tried but proved unsatisfactory (see 
also Sections 6.4 and 9.2). 
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Fig. 19 — Graphical summaiy of the HTB fit, (a) curves of y' vs L for constant 
x, (b) curves of y' vs x for constant L, (c) contours of constant y' in x,L space. 
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also small compared to the scatter in Fig. 10. This implies that a 
substantial fraction of the scatter may be associated with the high 
correlation between S and x c on the L-slice fits. Further consideration 
of standard errors and correlations of the fitted coefficients and 
detailed statistical evaluation of the fit is deferred to Section VIII. 

Fig. 19 presents a graphical summary of the function y'(x, L). 
Part (a) of the figure shows y' vs L for (several) constant x. Physi- 
cally, these curves correspond to values of the intensity of radiation 
vs L for constant magnetic dipole latitude, because x — constant 
implies A = constant. The nesting of the curves in Fig. 19(a) is a 
consequence of the fact that G'(x; x c , S) decreases monotonically 
with x [see (2) and Fig. 19(b)]. The shape of the curves changes 
smoothly with L, and the position of the maximum shifts smoothly 
toward higher L as the value of x (and therefore A) increases. 

The nesting property docs not hold for plots of y' vs x at constant 
L. This general consequence of the existence of a maximum in A'{L) 
is displayed in Fig. 19(b). All the curves in Fig. 19(b) have similar 
dependences on x. 




Fig. 19 — (continued) 
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Y (OBSERVED 



Fig. 20 — The value of y' computed from the HTB coefficients of Model I vs 
the observed value, Y, for the 960-point HTB sample. 

Fig. 19(c) contains contours of constant y' plotted in x,L space 
and completes the graphical summary. The contours surround the 
point x = 0, L = 1.46 at which the peak intensity occurs. 

7.3 Evaluation of Fit to the HTB Sample 

A summary indication of the quality of the fit of the 9-coefficient 
Model I to the HTB sample is given in Fig. 20, in which the fitted 
(computed) value, ?/, is plotted against the corresponding observed 
value, Y. The solid straight line would represent the case of a perfect 
fit. This is impossible on the basis of a model using only x,L 
coordinates since different Y values were observed for the same x,L 
pairs. It is seen, however, that the scatter of the plotted points about 
the line of perfect fit is reasonably uniform and that the horizontal 
width of the "scatterband" is roughly constant over the entire range of y'. 
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In the following subsections, the quality of fit to the entire body of 
HTB data is scrutinized, using many of the procedures used in the 
previous section to evaluate the CB fit. 

7.4 Evaluation of Fit on Equator 

The HTB fit along the equatorial boundary is displayed in Fig. 21. 
The points are the values of observed Y plotted against L for all HTB 
data for which ^ x < 0.037 (i.e., X < 1°), and the plotted curve is 
A'(L), defined in (6), using the HTB coefficients of Table IV. Comparing 
Fig. 21 with Fig. 11, it is seen that most of the excess scatter has been 
eliminated. The curve in Fig. 21 does not deviate noticeably from the 
center line of the points (except for 1.5 < L < 1.6, where the curve is a 
trifle high and for L « 1.95, where the curve is a trifle low). 




Fig. 21 — All the HTB data for x < 0.037 (i.e., within 1° of the magnetic in- 
variant, equator) and the equatorial value estimated from the HTB coefficients 
plotted against L. A' and Y are in units of ( counts/sec) J/a . 
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In Fig. 8 the solid curve, which is A'{L) calculated from the HTB 
coefficients, may be compared with the dashed curve, which is A'{L) 
calculated from the CB coefficients. The HTB fit gives higher 
equatorial values for y' when L is less than ~1.9, as might be ex- 
pected from the fact, displayed in Figs. 16(a) and (b) and discussed 
in Section 6.8, that the HTB data select the higher values of Y for 
1.4 < L < 1.6. For L greater than ~1.9, the equatorial values of the 
HTB fit are somewhat lower than those of the CB fit; however, there 
is no equatorial data for L > 1.95, and the comparison of the fits is 
not meaningful in this region. The points in Fig. 8 are estimates based 
on CB, not HTB, data and are not immediately pertinent to the solid 
curve. 

An estimate of the standard error of the fitted equatorial function 
A'{L), based on the HTB sample, is plotted as a function of L in 
Fig. 22(a) (see Section VIII for details). The standard error of 
A'{L) is typically less than one percent in the range of L (1.15 < 
L < 1.95) over which equatorial data are available. Error bars of 
this size would hardly be visible in Fig. 21. For the same values of L, 
the standard errors of A'(L) derived from the HTB fit are sub- 
stantially smaller than those from the L-slice fits listed in Table II. 
As might be anticipated, the percent standard error of A'{L) in- 
creases as the minimum x values of available data increases with 
increasing L beyond L = 2. This increase to a value of 10 percent at 
L = 3 reflects increasing uncertainty in the extrapolation of the fit. 
Note that the curves in Fig. 8, which represent the equatorial values 
of CB and HTB fits, differ, in general, by substantially more than two 
standard errors and the difference is certainly "statistically signi- 
ficant." 

7.5 Evaluation of Fit at Cutoff 

Figs. 23(b) and (d) show the positions, in x,L and R,\ coordinates, 
at which zero counts were observed during an 11-second counting 
interval. Figs. 23(a) and (c) are corresponding plots for one count 
(one, two, or three counts for L < 1.5) per counting interval. Only 
HTB data are plotted, and the density of points at high L has been 
reduced to improve the clarity of the display. 

Judgments regarding the quality of the fit are made, once again, 
with reference to the well-defined band of one count, rather than in 
terms of the more nebulous c utoff. The soli d lines in Figs. 23(a) and (c) 
are the loci of y'(x, L) = Vl count/11 sec, using the HTB coefficients 
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Fig. 22 — The standard deviation of A, a A , and the standard deviation of x c , 
a Xc , as functions of L. Units of a A and a Ze are the same as the units of A and x c , 
respectively, (a) Model I. (b) Model II. 
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in the model. These lines represent the data well. All hough the fit 
appears uniformly good in the x,L representation, a slight weakness 
near the "corner" at X ^ 40° is displayed sensitively in the R,\ plot 
(see also Fig. 12). 

The dashed lines in Figs. 23(a) and (c) show the locus of the fitted 
cutoff function, x c (L), calculated from the HTB coefficients. Error 
bars indicating excursions of one standard error in x c (L) are shown at 
two places on Figs. 23(a) and (c). The standard deviation of x e (L) 
as a function of L has been estimated (see Section VIII), and is plotted 
in Fig. 22(a). This standard error is smaller than those produced by 
the L-slice fits at corresponding values of L (see Table II). 

The values of x r (L) for the HTB and CB coefficients are plotted in 
Fig. 9. Although there is no discernible difference between the two 
curves in the figure for L < 2, the difference between the tabulated 
values exceeds twice the standard error (which is very small) over 
much of the range of L. The two sets of coefficients thus lead to results 
which differ in a "statistically significant" manner. For L less than 
^2, the significance of the standard error is more readily understood 
when it is interpreted in terms of the altitude of the cutoff. This is 
done in Section XI. 

Beyond L & 2, the values of x c for the CB and HTB coefficients 
diverge noticeably, compare Figs. 12(a) and (c) with Figs. 23(a) and 
(c), respectively. The magnitude of this divergence is quite sensitive 
to the method used in selecting the samples to be fitted. As has been 
discussed, the concept of a cutoff is not well defined in the context of 
these measurements for L > 2. The uncertainty is reflected in the 
rapid rise in the value of the standard error of X C (L) (see Fig. 22(a)] 
as L approaches 3. The significance of this rise may be more readily 
appreciated by referring once more to the error bars associated with 
x c (L) in Figs. 23(a) and (c). 

The partitioning of the data on the basis of electrical bias and tem- 
perature, and the procedure chosen for selecting the sample to the 
fitted, introduce statistically significant differences between the values 
of x e {L) obtained from the HTB and CB fits, as well as the more 
readily anticipated significant differences in the values of A'(L). 

7.6 Standard Error of Fitted Value 

The standard error for y'(x, L) is relatively constant, ranging be- 
tween 0.01 and 0.04, except close to x c (L). It should be understood 
that this standard error is based on the fit to the HTB sample, and 



1366 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1907 




Fig. 23 — All positions for the HTB data in R, X space (a) and x, L space (c) 
at which one count (one, two, and three counts for L < 1.5) was observed in an 
11-second counting interval, and all positions in It, X space (b) and x.L space 
(d) at which zero counts were observed in an 11-sccond counting interval. The 
solid lines are the loci of positions at which the HTB coefficients estimate one 
count in 11 seconds. The trace R = 2.0 R,, which explains the absence of data 
in the lower right-hand corner of the x, L plots, appears in part (d). The dashed 
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Fig. 23 — (continued) 

lines are the loci of the cutoff function x,(L) or R C (L) calculated from the 
HTB coefficients. The cluster of points near R = 1.1 R. and X = 20° in part 
(b) of the figure is data acquired by the telemetry station at Woomera, Aus- 
tralia. They represent observations made near perigee when the satellite was be- 
low the bottom edge of the proton belt, which is high over the western Pacific 
Ocean. 
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thus applies to the estimate of the average value of y and does not give 
the standard deviation of a single pred icted observation. The latter 
would be in the neighborhood of VOM = 0.2 (where 0.04 is approx- 
imately the MSE, see Table IV). 

Contours of constant percent standard error in the counting rate, 
y-, are shown by the curves in Fig. 24(a). For /, < 2 the standard 
error is less than 2 percent except close to the cutoff, where the value 
of y* is falling fast. (Near the cutoff, the standard error in x a is more 
informative.) In the absence of a fitted function, it would be neces- 
sary to average between about 30 and 300 observations to achieve a 
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Fig 24 — Contours of constant percent standard deviation in the counting rate, 
y\ calculated from the fits to the HTB sample and plotted in x,L space, (a) 
Model I. (b) Model II. 
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Fig. 25 — HTB residuals of F (i.e., F — i/ calculated from the HTB coefficients) 
plotted against L. The arrows indicate ± the approximate standard deviation if 
F 2 were Poisson distributed. 

standard deviation as small as 2 percent. As discussed in Appendix 
B.4, the estimates of the standard deviation based on the HTB sam- 
ple are conservative and (if there were no biases in the model) the 
values that apply to the 40,000 HTB points might he smaller than 
those in Fig. 24(a) by a factor as large as 6. 

The values in Fig. 24(a) are for relative counting rates (or fluxes) 
and do not include the uncertainty in the absolute calibration of the 
instrument noted at the end of Appendix A. Other discussion is given 
in Sections 9.4 and 12.2 and Appendix B.4. 

7.7 Behavior of the Fit on Several L-Slices 

Using the HTB coefficients, values of <M.r) were calculated for 
L m = 1.35, 1.805, 2.0215, and 1.79. The results are plotted as the 
heavy solid lines in Figs. 4 to 7. Recall that the points in these figures 
are not all HTB points. In general, the HTB points are those with the 
higher values of Y, although this may not be the case at L zz 2.2 
because of the temporal effects discussed in Section X. The four 
figures also allow further appreciation of the difference in results 
between CB fit and the HTB fit produced by the partitioning of the 
data and the refinement of the procedure by which the sample was 
selected. 

7.8 Residuals in x f L Spaee 

The residuals, Y — //, were computed for all the HTB data using 
the HTB coefficients. Fig. 25 is a plot of residuals against L, and 
Figs. 26 and 27 are plots of residuals against .r, in the indicated 
/./-ranges. These plots are analogous to Figs. 13 to 15, and as they 
display properties similar to the earlier figures, the discussion of 
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Fig. 26 — HTB residuals of Y (i.e., Y - y calculated from the HTB coefficients) 
plotted against x for 1.40 < L < 1.60. The arrows indicate ± the approximate 
standard deviation if F 2 were Poisson distributed. 

Section 6.6 applies. In particular, there is little indication of a de- 
pendence of the residuals on the magnetic coordinates. Moreover, the 
residuals in Figs. 25 to 27 are more closely clustered about zero than 
those in Figs. 13 to 15, confirming the fact that there is less scatter 
in the HTB data. This reduction in the scatter is especially marked 
in the neighborhood of the peak of the radiation belt (near x = 
between L = 1.4 and L = 1.6, Fig. 26) . 

7.9 Mean Square Residuals in x,L Space 

A breakdown of the mean square residuals (MSR) by i-ranges 
for the fit to the HTB data is given in Table III. This analysis is 
analogous to that presented in Section 6.7 for the CB fit. For the 
Group I data the MSR for the overall fit (1.1 < L < 3.0) is about 
(1.5) (0.023) = 0.036 and the largest entry under HTB Group I is 
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0.059. The anomalous trend of the MSR near L = 1.4 evidenced in 
the fit to the unrestricted data (see Section 6.7) has been largely 
eliminated. The overall MSR for the Group I data has been reduced 
by 15 percent. 

The breakdown of the MSR by L-ranges is not a particularly 
refined test of the quality of the fit. This index is based on essentially 
all the HTB data and, because the averaging procedure is blind to 
the distribution of data within L-ranges, favors results that fit best 
where the density of data is high. As the HTB sample was selected 
using criteria dependent on the preliminary fit to the data and does 
not necessarily favor x,L regions in which large quantities of data 
were acquired, the results of fitting this sample does not produce the 
lowest obtainable value of MSR for all of the HTB data. Examina- 
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Fig. 27 — HTB residuals of Y (i.e., Y - y calculated from the HTB coefficients) 
plotted against x for 1.S5 < L < 1.90. The arrows indicate ± the approximate 
standard deviation if Y 2 were Poisson distributed. 
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tion of the MSR in x,L cells shows the effect of the sample selection 
procedure on the MSR in L-ranges. Appendix C contains further in- 
formation and analysis of MSR in x,L cells. 

Model I with the HTB coefficients, provides a summary of the 
HTB data that, in the light of the many sources of variability and 
measurement errors, reasonably approaches the limit set by expected 
statistical fluctuations. 

7.10 Sources of Variability in the Data 

The residuals for the HTB data are now examined to see whether 
further identifiable sources of variability may be associated with 
them. Possible sources are: instrumental effects, errors in the ephemeris 
of the satellite, errors in the description of the magnetic field, telem- 
etry errors, fluctuations in the length of the counting interval, de- 
ficiencies in the model, and temporal variations. While all these must 
make some contribution to the MSR, the interrelationships among 
the coordinates discussed in Section 6.8 and the small size of the 
individual contributions, make positive identifications very difficult. 
We have not attempted to examine in detail the large number of 
small, apparently systematic, deviations discernible on the residual 
plots, although some of these may be "statistically significant." In- 
stead we have restricted our study to effects which are readily ap- 
parent on the residual plots. Where the observations are dense, an 
effect would be glaringly apparent if it introduced a shift of ^ 0.05 
in the local mean of the residuals. (This corresponds to a change of 
about 1.2 percent in flux at the peak of the proton intensity, and 
about 12 percent when the flux is a hundredth of its peak value.) 

Instrumental effects are associated with temperature, bias voltage, 
radiation damage, and imperfections in the omnidirectional char- 
acteristics of the detector. Restricting the range of temperature and 
bias voltage removed the major fraction of the instrument.il effects 
associated with these variables. Directional effects in the detector 
might show up when the residuals are plotted against y, the angle 
between the spin axis and the local magnetic field vector. However, 
no dependence was observed, indicating that the detector is effectively 
omnidirectional. Radiation damage, though technically an instru- 
mental effect, is more logically treated with temporal variations. 

Examination of plots of residuals versus various geographic co- 
ordinates did not reveal any systematic dependencies. In view of the 
small excess of the MSR over expectation for a random Poi.s.son 
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process, and the existence of other sources of error, it seems reason- 
able to conclude that the ephemerides were computed with sufficient 
accuracy for this analysis. 

The plots of residuals against the geographic coordinates as well 
as against x and L values were used to judge the quality of the co- 
efficients used to calculate the magnetic coordinates L and x. No 
systematic effects that can be attributed to flaws in the coefficients of 
the magnetic field were discerned. Nor is there any indication, in the 
form of excessive scatter of the residuals, that L is an imperfect 
coordinate in any part of the region of space covered by these data. 

Gross telemetry errors and those that occur in conjunction with 
noise bursts are easily identified and have been discarded. There 
remain telemetry errors that are indistinguishable from good data 
on a point-by-point basis, and these erroneous data must make some 
contribution to the scatter. As noted in Section 8.1, the distribution 
of the residuals has been looked into and they are found to be veiy 
well-behaved. However, it is not possible to make any quantitative 
estimates of the contribution of the remaining telemetry errors to the 
MSR. 

Temporal variations are an important source of variability, and 
Section X is devoted to their analysis. 

VIII. STATISTICAL CRITIQUE OF MODEL I. 

This section presents further information on statistical evaluation 
of the Model I fit. (Some background concerning relevant statistical 
techniques is given in Appendix B.) While confirming the very satis- 
factory performance of Model I in fitting the data, as presented in 
Section VII, some unsatisfactory aspects are uncovered and several 
defects of the model are pinpointed. The rectification of these defects 
is effected by use of Model II, discussed in Section IX. 

s.i Fit of Model I to Ike 900-point HTB Sample 

The analysis of variance for the fit of Model I to the 960-point HTB 
sample is shown in Table IV. This gives various partitionings of the 
total sum of squares (about 0) of the 960 observations (on the square 
root of counting rate scale). Table IV indicates the relevance of the 
model to the data in terms of its statistical effectiveness. Fitting the 
nine coefficients of the model accounts for more than 99.3 percent of 
the total sum of squares of the observations, leaving less than 0.7 per- 
cent associated with "error" or lack of fit. On a per degree-of-freedom- 



1374 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1967 

basis, the ratio of mean square for "fitted model" with 9 degrees of 
freedom to mean square for "error" is over 16,000. 

Of course, simply fitting the mean of all the data accounts for a 
sum of squares of 2121.2 of the total of 5374.7. Of the remaining "cor- 
rected" total sum of squares about the mean of 3253.6, the part of the 
model "orthogonal" to the mean accounts for 3218.9, i.e., approxi- 
mately 98.9 percent (so that the squared multiple correlation coefficient, 
R 2 , is 0.989). The corresponding ratio, mean square for the model with 
(9-1) =8 degrees of freedom to mean square for error, is over 11,000. 

It is worth emphasizing that the sample selection process which was 
used (see Section 7.1) is such that fitting the sample is, on a per ob- 
servation basis, a more challenging problem than it would be for the 
entire body of data (see Appendix B.3) . 

A summary graphical indication of the appropriateness of the fit is 
given in Fig. 20 which shows the fitted value plotted against the ob- 
served value. A perfect fit (essentially impossible here with any model 
based on x,L coordinates because different integral values of Y are 
observed near the same x,L point) would be the diagonal straight line 
shown. Deviations from fit should be gauged as horizontal spread about 
the line, since the observed quantities are plotted as abscissa, and are 
seen to be reasonably uniform throughout. 

Incisive indication of the quality of fit was provided by various 
plots of residuals (against L, x, y, time, etc.). Some representative 
plots over all the HTB data are shown in Figs. 25 to 27 and Figs. 41 to 
43. 

As a further examination of the adequacy of the fit to the selected 
HTB data, normal and half-normal probability plots (see Appendix 
B.8) were prepared for the 745 residuals comprising the subset of the 
960-point HTB sample for which x < x c (L). These plots are shown in 
Figs. 28 and 29. 

Fig. 28 does display a generally good linear configuration indicating 
that the residuals may reasonably be regarded as a sample from a nor- 
mal distribution. There is no suggestion of general asymmetry or other 
distributional peculiarities. There are perhaps three values which are 
statistically "too large," but not wildly so. Indeed, the plot is remark- 
ably well-behaved and reassuring. 

From some points of view, it is useful to consider the. statistical be- 
havior of the residuals without regard to their sign. Fig. 29 is a plot 
of the ordered absolute residuals against standard half-normal (folded 
standard normal) quantiles. This presentation is more focussed and 
sensitive to a statistical overabundance of large absolute residuals. The 
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plot is also very well-behaved, with indication of the same three overly 
large values. 

The reason for omitting from these plots all residuals from points 
for which x > x (L) is that, for those, the predicted value y is and, 
in the great majority, the observed Y was 0; hence, the residual is 0. 
Since it was exactly this information which determined the estimate 
x c {L) and since one could hardly expect a collection which includes 
about Ys zeros to behave like a normal sample, these points were omit- 
ted 

From either Figs. 28 or 29 one can estimate a slope of about 0.21, 
which is an estimate of the standard deviation of the (counting rate)* 
observations, clear of the confounding influence of the nonvariance- 
stabilized very low counting rate observations, since observations for 
x > x c (L) have been omitted. The corresponding variance estimate, 
0.044, clearly exceeds that from the Poisson approximation, 0.023, 
and also is greater than the pooled value for the AISD(y), 0.039, 
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Fig. 28 — Normal probability plot of residuals from fit of the model to the 
HTB sample. 
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(Appendix C) the overall HTB data MSR(F), 0.038, (Appendix C) 
as well as the MSE(F) from the fit to the 960 points, 0.036, (Table IV). 
This is as one would expect, since the variance estimate from the slope 
of Figs. 28 and 29 is not downward biased by the zero (and Vl/11) 
residuals from the very low counting rate observations for x > x r (L), 
while the other quantities are so biased. 

The excess of the variance estimate of 0.044 over the Poisson value 
of 0.023 may be due to any or all of several factors, including: (t) the 
non correctness of the Poisson assumption, (it) temporal variations in 
the radiation belts or the detection equipment, (tit) measurement 
errors or computational biases in time record, ephemeris or magnetic 
coordinates, etc. (to) noise bursts— the outlandish values were detected 
and discarded, but the general effect must be an upward bias on varia- 
tion, and (v) inadequacies in the model, including analytic form and 
coordinates employed. 
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Fig. 29 — Half-rioimal probability plot of absolute residuals from fit of the 
model to the HTB sample. 
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8.2 Statistical Measures Over All the HTB Data. 

An extensive presentation and comparison of various functions of 
the residuals over all the HTB data is given in Appendix C. Those re- 
sults provide (i) an empirical justification for the use of the square 
root transformation; (j?l a strong indication that the fit attained by 
Model I cannot he improved very much in the least squares sense over 
all the HTB data; {Hi) information on the extent of "unevenness" of 
the cell-construction process by which the 960-point HTB sample was 
selected; and (iv) some indication of differential effectiveness of fit of 
Model I to the data for different x,L regions. 

8.3 Statistical Properties of Estimates oj the Coefficients and Coefficient 
Functions. 

The least squares estimates of the nine coefficients of Model I fitted to 
the 960-point HTB sample are given in Table IV, with their approxi- 
mate standard errors and pairwisc correlations.* These provide the 
information needed to obtain estimates and standard errors for func- 
tions of the coefficients; e.g., y'i.r.L) , or A'{L) , or the value of the max- 
imum counting rate, or the position in space at which the intensity of 
high energy protons is maximum, etc. (See Appendix B for the neces- 
sary formulae.) 

Some of the pairwise correlations in Table IV are exceedingly high. 
This may be due, in general, either to an unfortunate "design" (i.e., 
the array of positions of observations in x,L space in this application) 
or to some inherent "coefficient redundancy" in the model, or to both 
such blemishes. Occurrence of such near-singularities can lead to prac- 
tical difficulty with the iterative fitting computation and/or make the 
individual coefficient estimates poorly determined. 

In the present model, only the coefficient L„ has a direct physical 
interpretation. Its estimate has a very small standard error and an 
entirely bearable correlation with the remaining coefficient estimates 
(all values of |«| < 0.5). Otherwise, physical interest centers mainly 
on the coefficient functions A'(L), x,.(L), and y'(x,L) whose estima- 
tion is considered in Sections 7.2, 7.4, 7.5, 7.6, and 8.4. 

For a given model and specified coefficient values, the matrix of ap- 
proximate correlations depends only on the array of data positions in 
x,L space. Thus, to check on whether the correlational problems might 


* A rescaling of the values of p. namely as the quantity a denned and moti- 
vated in Appendix B.5. is also given in Table IV. The coefficient of dependence 
a has more nearly (he behavior of a "linear utility function." 
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be due to inadequacy of the practically available (selected) array, a 
correlation matrix was computed using an 'ideal' x,L array, namely 
the 1034 values of (x,L) corresponding to the division of x,L space 
described in Section 7.1 and Appendix B.3. While some minor improve- 
ments in some of the correlations were noted, the changes were small. 
Thus, it would appear that the main reason for the high correlations 
is in fact some "coefficient redundancy" in the model. 

Inspection of Table IV indicates that the very large correlations are 
associated with some of the parameters of the A'{L) function, namely 
flit «2, as, and 77 for all pairs of which |p| > 0.99 (i.e., \a\ > 0.90). 
Moreover, it will be seen in Section 8.5 below, that the present param- 
eterization of the model leads to a markedly large indication of non- 
linearity and there is reason for believing that this is largely due to 
the same subset of coefficients. The combination of both defects stimu- 
lated development of Model II which overcame them (see Section IX). 

8.4 Estimates of Functions of the Coefficients 

The estimates of the coefficient functions A'(L) and x c {L) have been 
discussed in Sections 7.4 and 7.5 and summarized in Figs. 10 and 11. 
Their estimated standard deviations, on a "pointwise" basis, are 
graphed in Fig. 22(a), while the approximate correlations of the esti- 
mates of A'(L), x c (L), and S, as functions of L, are shown in Fig. 30(a). 

Despite the near-singularities (i.e., \ p \ near 1) in the estimates 
of some of the individual coefficients of A'{L), it is seen that the estimate 
of the square root of the equatorial counting rate provided by A'(L) is 
well-determined over the entire L range. The standard error varies 
between approximate limits of 0.018 and 0.040, nonmonotonically, and 
these values are typically less, sometimes by a factor of 5 or more, 
than the standard errors from the corresponding Z,-slice estimates 
(see Table II) reflecting in part the statistical gain from the simul- 
taneous two-dimensional fit. 

For x c (L), the standard error is less than 1 percent over much of 
the range of L, rising to 3 percent for large L values where the data 
are statistically inadequate. 

The three correlation functions p A , Xc (L), p A , s (L), and p SiXe (L), for 
the estimated coefficient functions A'(L), x c (L), and S, are plotted in 
Fig. 30(a) (see Appendix B.4 for formulae). In general, these correla- 
tions are small (\ p \ < 0.5, | a | < 0.12). The statement applies to 
the correlations involving A'{L) despite the very high correlations among 
individual coefficients. The generally low correlation between A'(L) 
and x c (L) is as intuitively expected since A'{L) is influenced mainly 
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Fig. 30 — Correlation coefficients of A with S, S with x , and A with x e , cal- 
culated from the fits to the HTB sample and plotted as functions of L. (a) Model 
1. (b) Model II. 
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by observations at small x while x c (L) is determined mainly by those 
at large x. The exception is near L = L , where p AlXc (L ) approaches 
1 as a result of the fact that the coefficient L is common to both func- 
tions and that the forms of A'{L) and x e (L) [see (4), (5), and (6)] 
require that both functions be zero when L = L . 

The statistical correlation between the fitted A and x c for the L-slice 
fits was always positive (see Table II), which is not the case for p A , Xt (L). 
This change in sign gives some indication of basic differences in be- 
havior between the results of the two-dimensional fit and the outcome 
of the collection of one-dimensional L-slice fits. 

The (A, S) and (S, x e ) correlations have the same signs in all cases. 
The magnitude of the correlations among A, x r , and S is larger for 
the L-slice fits (see Table II) than for the HTB fit at corresponding 
values of L [see Fig. 30(a)]. This is very noticeable for L greater than 
»1.7. particularly for the large correlation between S and x c . It is 
these large correlations which make it difficult to obtain reliable L-slice 
estimates of x c or S when L„, > 2 (see Fig. 6) or when the distribution 
of the data within an L-slice is poor (see Fig. 7). 

8.5 Nonlinearitij Indices and Dependence of Estimates 

Appendix B.5 discusses the use of the sum of squares function (i.e., 
sum of squares of differences between observed value and "fitted" 
value, as a function of proposed coefficients) as an indicator of the 
joint dependence and behaviour of the coefficient estimates and the 
fact that the extent to which the contours of the sum of squares func- 
tion arc approximated by a certain family of ellipsoids provides a meas- 
sure of linearity of the model. 

Fig. 31 shows 4 of the 36 pairwise projections of the 9-dimensional 
ellipsoid, whose size would correspond to a "0.99 joint confidence co- 
efficient" as discussed in Appendix B.5. The axes are scaled in each 
case according to the standard error of the coefficient. The orientation 
and shape of the ellipse corresponds directly to the sign and magnitude 
of the correlation, P , or its transform, «, for the pair of coefficients. 
Thus, for example, Fig. 31(a) shows the projection onto the ai-a a 
plane. The resulting very narrow positively inclined ellipse corresponds 
to a very high positive correlation of «i, Og (p = 0.9995, « = 0.97). 
(The 45° inclination of the graphed ellipses is a result of scaling the 
axes by their standard errors.) Part (b) of the figure shows a narrow 
negatively inclined ellipse for the case of rather large negative correla- 
tion between a 3 and v estimates. Parts (c) and (d) illustrate results for 
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small and negligible correlations between L , r 3 and r 8) S, respectively. 

At various positions on these ellipses there appear numbers which 
are ratios of the actual sum of squares at that "point" to the minimum 
sum of squares. The computation of the actual sum of squares is done 
for the coefficient values corresponding to the point on the 9-dimen- 
sional ellipsoid which projects into the point on the plotted ellipse. 

If, in fact, the coefficients occurred linearly, all of these numbers on 
all of the pairwise ellipses would be constant and in the present case 
would have the value 1.023 corresponding to a sum of squares of resid- 
uals of about 35.47. As a basis for judging the actual values and their 
variability, the following table gives values which this ratio would 
have, if the coefficients did occur linearly, for various joint (9-dimen- 
sional'l "confidence coefficients:" 



Con]. Cot 


\ff. 


C 


ontour Ratio 


0.90 






1.015 


0.95 






LOIS 


0.99 






1.023 


0.999 






1.029 



In view of the variability of the actual ratios in Fig. 31, and of the 
extent to which some depart from the values in the above table, it is 
clear that in the present form of the model the coefficients behave 
jointly in a markedly nonlinear fashion even in a relatively small 
neighborhood around the least squares estimate. 

Inspection of the entire set of (9) (8)/2 = 36 pairwise plots strongly 
suggests that a major part of this nonlinear behavior derives from the 
coefficients a,, a- 2 , a a , and -q of the A'(L) part of the model. These also 
arc the coefficients whose estimates exhibit the undesirably high cor- 
relations which have been shown above to be due mainly to a "coeffi- 
cient redundancy" in the model. 

Direct interpretation of the ellipses in Fig. 31, as indicating inter- 
dependence of the coefficient estimates, depends heavily on the appro- 
priateness of the linear approximation in the neighborhood of the least 
squares estimate. Since the nonlinearity index is in fact distressingly 
large one must be cautious in interpreting the ellipses or their asso- 
ciated correlation or dependence coefficients. 

8.G Summary Statistical Criticisms of Model I. 

Model I, with coefficients determined by fitting to the 960-point 
HTB sample, has been shown to provide a very good fit both to the 
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Fig. 31 — Examples of projections of the approximate "0.99 joint confidence 
region" for the estimates of Model I. (Axes are scaled by standard errors.) 
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sample and to the entire body of some 41,000 HTB observations. 
Moreover, the interesting coefficient functions y'(x,L) , A'(L) , and x c {L) 
have stable statistical properties as has the physically interpretable 
coefficient L Q . 

However, the model has two statistical defects: Firstly, although 
the model gives an extremely good fit to the data, the parameters 
fli , (h, Oa, and 77 of the A'(L) part of the model have exceedingly high 
mutual correlations (see Table IV), and these were shown not to be 
due to an obviously defective design. Secondly, the model coefficients 
exhibited distressingly high nonlinearity of behavior even within 
rather close neighborhoods of their least squares estimates, with 
grounds to suspect that this was caused by the ai, 03, &&, y group of 
coefficients. In addition, most of the coefficients of Model I do not 
have any directly meaningful physical interpretation. 

The modifications which led to Model II, as discussed in the 
following Section IX, overcome these defects of Model I while re- 
taining all its virtues. 

IX. THE MODEL II FIT TO THE HTB DATA 

This section presents the statistical analysis of the HTB data 
using Model II, a modified version of Model I. The emphasis in the 
presentation is on comparisons of Models I and II. Since it is shown 
how very closely the fit of Model II approximates that of Model I, 
such aspects as the direct presentation of Model II residuals overall 
the data are unnecessary, and hence omittted. 

9.1 Model II 

The definition of Model II has been given in Section 4.6, together 
with a discussion of the physical interpretation of its coefficients and 
its mathematical relation to Model I. Specifically, the 8-coefficient 
Model II constitutes a specialization and reparameterization of the 
9-coefficient Model I. Thus, it follows that the minimum sum of 
squares in fitting Model II to any body of data can not be less than 
that from fitting Model I, though this may not be true of the mean 
square error. 

The evolution of Model II from Model I did not arise from any 
simply described systematic process, as is indeed true in other aspects 
of this study. Once the basic achievements of Model I were estab- 
lished it was then opportune to focus on major remaining defects. The 
character of these defects strongly urged elimination of one or more 
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coefficients in conjunction with a nonlinear reparameterization of 
the coefficients. The solution achieved was arrived at by empiricism, 
persistence and good luck. 

The remainder of this section documents the assertion that Model 
II retains all the virtues of Model I while overcoming its defects. 

9.2 The Fit of Model II to the 960-point HTB Sample. 

The analysis of variance from fitting the 960-point HTB sample 
by means of the 8-coefficient Model II is given in Table V. As ex- 
pected, the residual sum of squares, 34.7126, of Table V exceeds that of 
Table IV, namely 34.6675. This difference is associated with the 
one-degree-of-freedom nonlinear constraint defined in (13). Thus, 
we see. that the sum of squares associated with the one-degree-of- 
freedom non-linear constraint is (34.7126-34.6675) = 0.0451 and this 
gives a ratio of less than 1.24 in relation to the mean square error 
of 0.03645. The value 1.24 corresponds to the upper tail 27 percent 
point of the chi-squarcd-with-one-degree-of-freedom distribution. 
The proportionate increases in the sum of squares for error is about 
0.13 percent and the increase in the mean square error is less than 
one part in 3000. Multiple R- = 0.989 is effectively unchanged. 

For the models of both Tables IV and V, the coefficient S is treated 
as constant with L. If Model II is modified so that S(L) = s n + SiL, 
then, fitting this 9-parameter version of Model II yields a sum of 
squares for error of 34.520. Thus, we would have a sum of squares of 
(34.713-34.520) = 0.193 associated with the "hypothesis" x, = 0. 
The main point of quoting this result is to indicate that these minor 
differences in the sums of squares for error are judged as unimportant 
in this context, even if under some highly formalized assumptions the 
distinctions are "statistically significant." 

Of greater interest and sensitivity are the following considerations: 
(0 the behavior of the residuals from Model II as functions of x,L 
and y; (it) the behavior of the differences between Models I and II; 
(iu) comparisons of the estimates of A'(L) of Model I and A"(L) 
of Model II | see (6) and (11) |; [iv) comparisons of the estimates of 
x,.(L) from the two models; (v) the pattern of correlations of the 
estimates of the eight Model II coefficients; and (vi) the indices of 
nonlinearity for the coefficients of Model II. 

9.3 Residuals of Model II Fit and Differences Between Models I and II. 
Figs. 32, 33, and 34 are plots of the residuals of the 960-point HTB 

sample from the fitted values of Model II against L, x and Y, re- 



1386 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1967 



W 

3 

Oh 
§ 
< 
OS 

a 



3c 

o 

o 
o 

p 



fa 



w 

a 
< 



c £ 


667.5024 
0.0365 


o v 

11 


5374.7321 

5340.0195 

34.7126 


■d 


OGOCN 
CO "O 
OS OS 


o 
** 




Total 

Model 

Error 



■s: 


ION 
N -+ 
(N O 

CO O 

do 


t- 


00 CO 

coos 
c©-< 
coo 

do 


t. 


CM -^ 

CO-* 
0<N 
•CO 

do" 

1 


C 


coco 
mco 

COO 
(NO 

do 


p 


MN 
NO 

io© 


4 


32 

COO 

Tt<0 

i-iO 


4 


COOS 
OSO 

<N o 

HO 

T-id 


4 


CNN 

CO* 
NCO 
OO 

cod 




-a 
— 

a 

-U 
to 

<B 

|E 





BQ 


— < O -H O -H O O 

(NOOOOOO 

ooooooo 

1 1 1 




v- 


iO 

NfflHHOlH CO 

OOOOBh O 

oooooo o 

1 1 1 1 1 1 




t. 


0.01 

0.19 

0.00 

-0.00 

-0.62 

-0.956 
0.098 




c 


-0.02 

-0.41 
0.00 
0.02 

-0.923 

0.793 

-0.168 


8 


E- 


0.02 

-0.11 

0.39 

0.214 
-0.021 
-0.107 
-0.010 




^ 


O iO OS 5f OS iO OS 
O O NOOHO 

do ©odd© 
1 1 




4 


■* i-H CO OS <N »C 

H OWOCONH 
O COtTCOiO-* o 

o oooooo 

III 1 




a 
^ 


O CO CO -H IN iH OS 
■* -ROCOCO t^O 

OOOOOOO 
1 1 






ft o ft 



protox'data from telstar i 1387 

spectively. These plots show no systematic structure and are quite 
similar to analogous plots for Model I. Furthermore, Fig. 35, showing 
the observed Y versus fitted y" for Model II, is as well-behaved as 
the corresponding Fig. 20 for Model I. 

Figs. 36, 37, and 38 show the deviations between the fitted Models 
I and II plotted against L, x, and Y, respectively. Of course these 
figures show a systematic structure since one is plotting the difference 
of two smooth functions. However, the actual differences are totally 
insignificant in the light of the data. (Note that the scale for Figs. 
36, 37, and 38 differs from that of Figs. 32, 33, and 34 by a factor of 
10.) 

Thus, on the basis of one less coefficient, Model II fits the data 
essentially as well as Model I, to which indeed it is a very excellent 
approximation. It has the merit that the physically arbitrary coef- 




Fig. 32 — Residuals (Y - y) from the fit of Model II to the 960-point HTB 
sample vs L. 
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Fig. 33 — Residuals (Y — y) from the fit of Model II to the 960-point HTB 
sample vs. x. 

ficicnts a, , a 2 , and a 3 of Model I have been replaced by A p and L p which 
do have direct physical interpretations. As will be detailed in the 
next subsection, Model II also has additional attractive statistical 
attributes. 

9.4 Coefficient Estimates 

Table V gives the least squares estimates of the eight coefficients of 
Model II together with their approximate standard errors, correla- 
tions and a values. The estimates are seen to be extremely well- 
determined. In particular, for the physically meaningful quantities 
A p , L , and L p the standard errors are about 0.4, 0.1, and 0.15 percent, 
respectively, while for the shape coefficients v and S they are about 
1 and 1.5 percent, respectively. 

Comparison with Table IV shows that the standard error has de- 
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creased for every coefficient which is common to the models. The 
most dramatic change is for -q for which the standard error diminished 
by a factor of about 8. 

The estimates of A'(L) and A"{L) are in very close correspondence 
as implied by Fig. 36. The comparison of Fig. 22(b) with Fig. 22(a) 
indicates that the standard error of A"(L) is uniformly lower than 
(but in general agreement with) that of A'(L). 

Entirely similar remarks apply to comparison of estimates of x e (L) 
from Models I and II, as also documented by Figs. 22(a) and 22(b). 

It has already been shown that the fitted values of y'(x, L) and 
y"(x, L) are in very close agreement. The pattern of contours of the 
percent standard errors of |.i/"(.r. L)] 2 , in Fig. 24(1)), shows that the 
standard error is everywhere smaller than the corresponding results 
for Model I, in Fig. 24(a). 
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Fig. 34 — Residual* ()' — y) from the fit of Model II to the 960-point HTB 
lample vs Y. 
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One of the most dramatic changes between Models I and II is 
indicated by comparison of the correlations in Tables IV and V. The 
very large correlations (\ P \ > 0.99, |«| > 0.9) among the A'(L) coef- 
ficients of Model I do not occur for Model II. Only the (r u r a ) and 
(r 2 , r 8 ) coefficient pairs of Model II have \a\ values above 0.5. This is 
inconsequential since these are physically arbitrary coefficients of a 
cubic polynomial. 

The correlations of A"(L),x (L), and S from Model II remain much 
like the corresponding results for Model I, as shown in Fig. 30. 

9.5 Nonlinearity Indices 

The further virtuosity of Model II is indicated by the behavior of 
the nonlinearity index shown for the examples of "confidence regions" 




Y (observed) 

Fig. 35 — The value of y" computed from the fit of Model II vs the observed 
value, Y, for the 960-point HTB sample. 
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Fig. 36 — Deviations between the Model-I fit, y', and the Model-II fit, y", vs 
L, for the 960-point HTB sample. 

in Fig. 39. (See Appendix B for general discussion and definition.) 
Specifically, it is seen that the numbers on the ellipses vary very 
little and this is true for all 28 of these ellipses. These numbers would 
be constant and all equal to 1.023 if the model were linear in the 
fitted coefficients. Comparatively, Model II does indeed behave in a 
reassuringly linear fashion. For sharp contrast, we may compare Fig. 
39 with Fig. 31, for Model I, in which the values range up to 1000 
around the 9-dimensional ellipsoid. 

The nonlinear behavior of Model I in relation to the linear be- 
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Fig. 37 — Deviations between the Model-I fit, y', and the Model-II fit, y", vs 
x, for the 960-point HTB sample. 
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Fig. 38 — Deviations between the Model-I fit, y', and the Model-II fit, y", vs 
}', for the 960-point HTB sample. 

havior of its specialized reparameterized version, Model II, is in- 
dicative of the reason for the high nonlinearity indices for Model I. 
Effectively, a p-coefficient model defines a constraining "surface" 
of p dimensions (p is 9 and 8 for Models I and II, respectively) in 
the rj-dimensional space of the observations (n is 960 in the present 
case). In a small neighborhood of the least squares estimate, this 
p-dimensional surface may or may not be planar. If the latter, one 
will obtain high indices of nonlinearity. If the former, then one will 
or will not obtain high nonlinearity indices according to whether 
the individual coefficient coordinates within the p-dimensional surface 
are or are not linearly behaved. 

It is likely that the 9-dimensional surface defined by Model I is 
indeed reasonably planar, but the coordinate system denned by the 
coefficients is highly nonlinear. 

The correlation and nonlinearity effects, it should be noted, are not 
in principle related. One can have very high correlations with linear 
models and very low correlations with very nonlinear ones. 

9.6 Summary Comments 

Model II has been presented and validated as an evolution of Model 
I. Though Model II represents the current recommended fit from 
this study, several aspects of its justification, and of other comparisons 
in this paper, are based on the Model I fit. For example, the statistical 
study of residuals over all the HTB data, discussed in various places 
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including Appendix C, is based on Model I. This hybrid attitude is 
entirely sound, since the range of deviation between Models I and II is 
small compared to the range of residuals from the fitted sample. 

Thus Model II provides a fit to the HTB in which the 8 estimated 
coefficients provide a "good description" of about 41,000 observations. 
The deviations of the fit from the data are within reasonable statistical 
fluctuations— variation in telemetered counting rates, orbital errors, 
observational errors, mapping-to-magnetic-coordinate uncertainties, 
etc. (See Appendix C.3). A number of the coefficients have physical 
interpretations and these are statistically well-determined and rela- 
tively uncorrelated. Model II, though nonlinear in the coefficients, 
behaves in a very linear fashion in the neighborhood of the least 
squares estimates. 

X. TEMPORAL VARIATIONS 

This section and the two to follow are devoted to discussion of 
some specific physical results of the analysis. 

Temporal variations are considered in three classes: diurnal (day- 
night), secular, and short term. Residual plots were used to study these 
effects. 

10.1 Diurnal Effects 

The HTB residuals were plotted against local time for various 
x,L regions. The HTB data are not well-distributed in local time 
near the magnetic dipole equator, making it difficult to draw firm con- 
clusions. However, no evidence of a diurnal variation was found. 

Specifically, to produce a change of about two percent in the 
average value of Y on the equator (x = 0) would require a diurnal 
shift in the radial position of the magnetic field line of about 0.01 
E,. at /, = 1.3f), and a shift of about 0.02 R r at L = 1.55, if there 
were no other effects. At these two positions, the value of y is large 
(// =r 8) and dij/dL is large, and a two-percent change in y would 
correspond to a shift in the mean of the residuals of ^0.16 between 
noon and midnight local time. An effect of this magnitude would be 
readily observable on the residual plots. 

Thus, it is unlikely that displacements larger than 70 km and 140 
km, at equatorial L's of 1.35 and 1.55, respectively, would escape de- 
tection, and these distances are offered as upper limits to the day- 
night changes of the magnetic field at the two positions. As both of 
these displacements are equivalent to a change in field strength of 
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Fig. 39 — Examples of projections of the approximate "0.99 joint confidence 
region" for the estimates of Model II. (Axes are scaled by standard errors.) 
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about 300 gamma (0.003 gauss), this particle experiment does not 
qualify as a sensitive indicator of adiabatic changes in the earth's 
magnetic field. 

10.2 Secular Effects 

The HTB residuals are plotted against elapsed time, in days, for 
1.85 < L < 1.90, in Fig. 40. It would appear that the average value of 
Y decreased between days 191 and 255. This decrease is exhibited in all 
parts of the belt where we have measurements during this interval. 
Between days 191 and 225, the orbit of the Telstar® 1 satellite did not 
take it into the central region of the belt {1.3 S L ^ 1.8, X ^ 10°}. In 
other regions the decrease in the average value of Y over this period is 
about ten percent. The extremes are two percent and 20 percent, but it 
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Fig. 40 — HTB residuals of Y (i.e., Y — y calculated from the HTB coefficients) 
plotted against time for 1.85 < L < 1.90. The arrows indicate ± the approximate 
standard deviation if I' 2 were Poisson distributed. 
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is not possible to separate out other variables which may be influencing 
the results. 

From the magnitude of this effect, it is clear that it must be con- 
tributing substantially to the MSR. A decrease of ten percent in the 
average value of Y corresponds to a decrease of about 20 percent in 
the flux. A fractional change in the flux which is independent of X 
and L cannot be distinguished from a change in the characteristics 
of the instrument. Among other possibilities, radiation damage or the 
decay of protons which might have been associated with the Star- 
fish high-altitude nuclear test of July 9, (day 190) 1962 might have 
produced the observed effects. Because of this ambiguity, we are 
unable to offer any well-founded interpretation of the time depend- 
ence of the data before day 225. For reasons to be noted shortly, 
ambiguities are also encountered when interpretation of the temporal 
behavior of data acquired after day 400 is attempted. In the inter- 
mediate period, the time dependence does vary with x and L. By 
using Fig. 40, which shows comparatively little fluctuation during 
this intermediate period, as a standard we are able to measure 
relative changes in the belt. The stretches of sparse data near days 
240 and 320 in Fig. 40 arc a result of the orbital configuration, there 
being less opportunity to acquire "high-temperature" data during 
these periods. The absence of HTB data between day 325 and 373 
was caused, as noted in Section 6.9, by the low bias condition that 
existed during that time. However, an examination of residuals from 
the CB fit between days 325 and 373 reveals nothing that vitiates the 
conclusions drawn from the HTB data in what follows. 

Residuals versus time-in-days have also been plotted for x,L cells 
of size 0.1 in L by 0.2 in r. Below L = 1.9 we find only one change 
with time within the sensitivity of our measurements, namely, a 
secular decrease between days 225 and 400 which occurs only near the 
ends of the field lines (x ^ x e — 0.2). We are unable to quantify this 
effect because, in order to see the droop above the noise, we need to 
collect residuals from a fairly sizable region of space. The term "sizable" 
means a region over which y changes so much that an average value of y 
in the region is not sufficiently representative to be used as a basis for 
computing a percent change in the flux. Fig. 41 gives an example of an 
x,L cell near the cutoff where this decrease may be seen. However, in 
the adjacent lower-a; region, Fig. 42, where the ability to discriminate 
absolute changes in the average value of Y is the same and the ability 
to discriminate percent change in the average value of Y is much greater 
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Fig. 41 — HTB residuals of Y (i.e., Y — y calculated from the HTB coefficients) 
plotted against time for 1.6 < L < 1.7 and 0.8 < x < 1.0. The arrows indicate 
± the approximate standard deviation if Y 2 were Poisson distributed. 

than for the region of Fig. 41, no corresponding secular decrease be- 
tween days 225 and 400 is evident. 

The droop in the residuals after day 400, which is noticeable in 
Fig. 42, is characteristic of many of the plots of residuals versus 
time-in-days. The widespread occurrence of this effect confuses in- 
strumental and "real" variations and introduces unresolvable am- 
biguities when attempts are made to identify the source of the droop. 

The observation of the general downward slope in Fig. 41 might 
be explained by a small decrease in x 0) which corresponds to a small 
increase in the altitude of the cutoff, between August 1962 and 
January 1963 on L-shells below 1.9." Alternatively, one might be 
observing the decay of the 55 MeV protons whose perturbation by 
the Starfish high-altitude nuclear test of July 10, (day 190) 1962 and 
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subsequent behavior have been measured by Filz 24 near the bottom 
of the trapped proton belt. There are too few data for us to attempt 
further interpretation of this qualitative observation concerning the 
secular behavior of x c . The number of points affected and the mag- 
nitude of the shift are too small for this effect to contribute interest- 
ingly to the MSR. 

10.3 Short-Term Effect 

The plots of the residuals versus time-in-days, for x,L regions, 
show a short-term fluctuation which is sufficiently singular to be re- 
ferred to as an event. This event is an increase in the average value 
of Y over the 30-day period which starts about day 280. It can be 
seen clearly in Fig. 43. The increase is discernible only for L > 1.9. 
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Fig. 42 — HTB residuals of Y (i.e., Y - y calculated from the HTB coefficients) 
plotted against time for 1.6 < L < 1.7 and O.G < x < O.S. The arrows indicate 
± the approximate standard deviation if }'- were Poisson distributed. 
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Table VI — Fractional Increase in Flux Between Days 280 

and 310, 1962. 



L^\x 


0.1 


0.3 


0.5 


0.7 


0.9 


<1.9 
1.95 
2.05 
2.15 


0.05 


0.07 


0.12 
0.37 
0.28 


0.20 
0.46 
0.33 


0.70 
0.90 



Table VI gives the fractional increase in the average counting rate 
(V 2 ) during this period at various values of x and L. By L = 2.25 
the change is barely observable and for L > 2.3 it has disappeared. 
The data acquired between days 325 and 373, which are not included 
among the HTB data because the bias voltage was low, were ex- 
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Fig. 43 — HTB residuals of Y (i.e., Y — y calculated from the HTB coefficients) 
plotted against time for 2.0 < L < 2.1 and 0.6 <x < 0.8. The arrows indicate 
± the approximate standard deviation if Y- were Poisson distributed. 
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amined; and there appears no reason to believe that there were any 
changes in the intensity of the >50 MeV protons for L > 1.9 during 
these 48 days. 

While it is not possible to be quite sure that we are observing a 
"true" temporal effect, it is difficult to contrive any alternate ex- 
planation. This event can be compared with the changes produced in 
the high energy proton distribution by the magnetic storm of Septem- 
ber 22, 1963, and observed with Relay l 25 and the Telstar® 2 satellite/' 
In both cases only L shells with values above 1.9 were affected, and the 
effect is more pronounced at higher .r's. However, the storm produced a 
decrease in flux whereas an increase was observed in 1962; the effects of 
the storm were more severe at larger L's, whereas in this event, a max- 
imum fractional change was observed near L = 2.05; and the effect 
of the storm was sudden, i.e., the flux decrease took place within 24 
hours, while the increase observed in 1962 was gradual and required 
a month to complete. Increases in flux having some of the features 
described here were observed with Explorer 7/ n However, it is dif- 
ficult to be certain that those increases were caused by protons with 
energies above 18 MeV, rather than electrons with energies greater 
than 1.1 MeV. 

The high-energy protons appear very stable over the seven months 
covered by our data. In particular, no effects associated with the 
USSR high-altitude nuclear tests of October 22, October 28, and No- 
vember 1, 1962, or the large magnetic storm of December 18, 1962 have 
been observed. 

In summary, changes through time in the observed values of the 
flux are generally less than 20 percent, although they may be larger 
in some regions of space. We have not been able to detect a diurnal 
effect. Often, secular changes are not separable from other variables, 
an exception being an apparent change in the position of the cutoff. 
An event which appears to comprise a measurable redistribution of 
the proton flux over an appreciable volume of space and period of 
time has been noted. We do not know whether the redistribution is 
in energy or space, and find no indication of the mechanism in the 
data. 

XL THE CUTOFF 

As discussed in Sections V, 6.3, and 7.5, the cutoff function, x,.(L) , is 
defined in terms of our instrument, model and fitting procedure. For 
L < 2, the value of x c {L) corresponds to the position on the given L 
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shell at which the omnidirectional flux is of the order of 1 proton/cm 2 
sec, more than three orders of magnitude below the highest flux in 
the belt. However, because the flux is falling so fast with x, this 
position is almost certainly very close to the place at which the flux 
becomes 0. The last statement is not true for L > 2. Here, although 
the value of x c (L) (the place at which y = 0) still corresponds to 
the point at which the limit of sensitivity of our instrument is 
reached, the position of x a is not so well-defined by the fit. In addition, 
one has only to examine Fig. 23 to realize that x c may be significantly 
removed from the value of x at which the flux falls to zero. 

The Model-I HTB coefficients of Table IV define the cutoff func- 
tion, and we have made use of a modification of R. H. Pennington's 
mirror trace program* to calculate the minimum altitude correspond- 
ing to x c (L) for L < 2.2. This inversion was accomplished using the 
Jensen and Cain magnitude field coefficients for I960, 13 the same set 
used to calculate x and L (see Table I). (Other sets of coefficients are 
available. 27 However, using the GSFC (7/65) coefficients 28 does not 
produce significantly different altitudes.) 

The minimum altitude is smallest in the Southern Hemisphere over 
the Atlantic Ocean. Fig. 44 shows the results in graphical form. The 
minimum altitude is ^270 km near the equator (L = L zz 1.13), 
decreases to a minimum of ^160 km at L = 1.6, and increases very 
rapidly thereafter. For L less than 1 .5, the standard error in altitude, 
derived from the standard error in x (see Fig. 22), is about 10 km, 
which is roughly the accuracy of the inversion procedure as we used 
it. The standard error in altitude for L > 1.5 is indicated by the 
dashed lines in Fig. 44. At L = 2, where the cutoff mechanism is only 
partially atmospheric, the standard error is nearly 50 km. 

The minimum near L - 1.6 in the altitude curve of Fig. 44 appears 
to reflect the existence of the South American magnetic anomaly. 
Although R C (L) [see (5)] increases monotonically with L for L > 1, 
the increase is apparently not fast enough to override the influence of 
the anomaly. This result is true for all the sets of coefficients pro- 
duced in many trial fits as well as for the HTB coefficients in Table 
IV. We have not yet carried out the obvious next step of averaging the 
atmospheric density over the orbital path of the protons to see 
whether or not the shape of Fig. 44 can be explained on the basis of 
present models of the atmosphere. 

Although the shape of the minimum altitude curve remains the 



* Kindly communicated to us by D. J. Williams. 
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.same, the value of the altitude is sensitive to the method of select- 
ing the sample (see Section 7.1). For example, the. minimum value of 
altitude calculated from the CB coefficients is 100 km (again at 
L = 1.6), 60 km lower than the 160 km calculated from the HTB 
coi'fncients. The weighting of the HTB sample emphasizes the high 
X data and gives better representation, and therefore a better ex- 
pectation of fitting well, near the cutoff. However, the Telstar® 1 satel- 
lite, with its eccentric orbit and relatively high (950 km) perigee, could 
not give detailed information about particles near the top of the atmos- 
phere, and this is reflected in the results of the analysis. 

In conclusion, the curve of Fig. 44 probably represents the quali- 




Fig. 44 — The minimum altitude reached by > 50 MeV protons as a function 
of L. This altitude is determined in geographic coordinates from the transform 
of x c (L). The dashed curves are ± one standard error. 
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tative behavior of the minimum altitude of the cutoff reasonably well, 
but the uncertainty in the value of the altitude is larger than a 
simple examination of the standard error plotted in the figure would 
lead one to believe. The implications of these results for the details 
of the cutoff mechanism have not been examined in detail; however, 
it is clear from the sudden upturn of the curve in Fig. 44 that the 
mechanism is principally atmospheric for L less than about 1.9 and 
principally nonatmospheric on higher L shells. 

XII. COMPARISON WITH OTHER WORK 

12.1 Introduction 

When making comparisons among the various high-energy proton 
measurements it is desirable that the results be extensive in time and 
space, reported in terms of omnidirectional fluxes at various positions, 
and that these positions be expressed in magnetic coordinates de- 
rivable from the B,L set. A list of some experiments which meet these 
desiderata is given in Table VII. 

Following a presentation of flux maps, comparisons among these 
experiments arc made with respect to the following features: the 
absolute intensity at one point in the belt, as close to the maximum of 
intensity as is practical; the intensity vs L in the equatorial plane; 
the behavior of the intensity on selected L shells; the flux near the 
top of the atmosphere, and the equatorial pitch angle distribution. 
Comparisons covering a larger range of proton energies have also 
been made by Vette'-" and Fillius. 20 

One of the difficulties encountered in making comparisons among 
the various bodies of data is that most of the results have been pub- 
lished in graphical form, rendering it necessaiy to scale numerical 
values from small plots, an inaccurate procedure at best. A welcome 
exception is the Explorer 15 data, which Mcllwain 18 has made avail- 
able by means of a series of interpolation functions in the form of a 
FORTRAN computer program. 

12.2 Tehtar 1 Flux Maps 

For this discussion, the Tehtar® 1 HTB results have been converted 
to omnidirectional flux, /, where ./ = 4-mf/g. (Note that the value of 
g derives from the assumptions of Appendix A regarding the energy 
spectrum.) This procedure provides an estimate of the flux of protons 
with energies between 50 and 130 MeV at positions, (.r, L), in mag- 
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MAGNETIC 
INVARIANT EQUATOR 0° 



Fig. 45 — Omnidirectional isoflux contours derived from the HTB coefficients 
and plotted in R,\ space. Dashes indicate extrapolation beyond the region in 
which data were acquired. Long dashes form contours of constant percent standard 
deviation. 



Label | A 1 B 


C 


D | E 




L 


Omnidirectional 5 X 10' 2 X 10 1 
flux (J) I 1 


1 X 10* 


5 X 10 2 2 X 10" 




1 X 10° protons/ 
cm 5 sec 



netic space on the basis of the presently provided model and fit to the 
HTB data. 

For ease of reference, Telstar® 1 HTB flux maps are presented in 
three, commonly used forms: Fig. 45 shows contours of constant flux 
in R,X coordinates; Fig. 46, contours of constant flux in B,L co- 
ordinates; and Fig. 47, log flux vs log B curves for various values of L. 
These three graphs give an overall picture of the particle distribution. 
In these figures, dashed lines are used to indicate the extrapolation of 
fitted values to regions not penetrated by the satellite. Note the way 
the geometry of the coordinate transformations affects the extrap- 
olated regions. In particular, the functional extrapolation in B,L 
coordinates gives much more curvature to the contours than might 
be anticipated. The difference between the functional and straight 
line extrapolation in B,L can be as large as a factor of 2 in the 
flux (a shift of 0.2 in L) at L = 3. Except for the region of the 
secondary local maximum in the flux near L = 2.2, this functional 
extrapolation compares surprisingly well with the measurements made 
on higher altitude satellites. ■"'■ 1S 

In the altitude range covered by the data, a single maximum is 
observed. This maximum in the omnidirectional flux of ^ 6 X 10 3 
protons/cm 2 sec is located on the magnetic equator at R = L = 1.46. 
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The intensity falls abruptly near the bottom of the belt (the top of the 
atmosphere) and decreases more gradually toward the sides and top 
of the belt. On a given L shell, the intensity is a maximum at the 
magnetic equator, and deceases monotonically as the distance from 
the equator increases. 

Neglecting the uncertainties in the calibration of the instrument 
(-25 to +50 percent), which are discussed in Appendix A and are 
mentioned in the next subsection, the estimated standard deviation of 
the estimate of J is less than 2 percent of J over much of the region of 
space discussed in this section. Smoothed contours of 1 percent, 2 
percent and 5 percent standard error are plotted as the dotted 
lines in Fig. 45. Near the cutoff, where the counting rate is falling to 
zero, the standard deviation in x (sec Figs. 44 and 22) is a useful 
indication of uncertainty in the flux. Other information concerning 



0.25 



3 0.15 
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Fig. 46 — Omnidirectional isoflux contours derived from the HTB coefficients 
and plotted in B,L space. Dashes indicate extrapolation beyond the region in 
which data were acquired. Labeling is given in Fig. 45. 
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Fig. 47 _ Omnidirectional flux on several L-shells derived from the HTB co- 
efficients and plotted against log B. Adjacent curves in parts (a) and (b) are 
slipped one decade in J. All curves rise from ./ = 1 proton/cm 9 sec. Dashes indi- 
cate extrapolation beyond the region in which data were acquired. 
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standard deviations may be found in Sections 7.4 to 7.6 and 8.4, Figs. 
22, 23 (a), 23(c), and 24. 

The equations denning Model II (see Section 4.6) and coefficients 
of Table V, together with the transformation equations among various 
magnetic coordinate systems, allow accurate relative flux values to 
be easily calculated in any coordinate system. 

12.3 Comparison of Absolute Intensities 

The solid curve in Fig. 48 is the fitted omnidirectional equatorial 
flux of 50-130 MeV protons measured by the Telstai-® 1 satellite. The 
points are fluxes observed on other satellites (Table VII) at the mag- 




Fig. 48 — Values of equatorial omnidirectional flux, for the satellites indicated 
in the legend, corrected to the energy range 50-130 MeV and plotted at the ap- 
propriate value of L. An integral power-law energy spectrum [see (17)] of ex- 
ponent —M, where M is given is a function of L by the dashed curve, was used 
in making the corrections. References are given in Table VII. 
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netic equator and corrected to 50-130 MeV by using a single-com- 
ponent integral energy spectrum of the form 

N(>E) « E~" . (17) 

The values of M at the magnetic equator are plotted as the dashed 
line in Fig. 48. These values were taken from Gabbe and Brown, 5 and 
are consistent with those of Brown, Gabbe, and Rosenzweig, 3 and also 
those of Fillius and Mcllwain, 34 and Freden et al, 3D where the data 
overlap. Because of uncertainties in the geometric factors of the de- 
tectors (see Appendix A) and changes in the belt with time (see 
Section X) , one might expect agreement only within a factor of about 
2. On this basis the agreement in absolute intensity is quite reason- 
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l«'ig. 49 — Values of equatorial omnidirectional flux, for the satellites and en- 
ergy ranges indicated in the legend, plotted against L. The dashed curve is the 
ratio of Telslar® 1 to Explorer 15 measurements. References are given in Table 
VII. 
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able. However, the Telstar® measurements are somewhat on the low 
side, and those of Imhof and Smith 32 {H 2 and H s on Fig. 48) are 
much higher than the other observations. 

The points represent measurements taken before, after, and during 
the Telstar® 1 experiment so it is unlikely that changes in the flux 
with time explain these differences. It is difficult to account for the 
discrepancies in absolute flux in terms of the spectral correction, 
unless more complex spectral forms than those of Appendix A are 
considered, because the comparisons are among results of detectors 
whose threshold energies are close to 50 MeV. The most likely sources 
of the differences are errors in absolute calibration. It follows that 
a good deal of caution should be exercised in drawing conclusions 
about temporal effects and energy spectra from measurements made 
with different instruments. 

12.4 Intensity vs L in the Equatorial Plane 

Fig. 49 is a plot of the omnidirectional equatorial flux for each of 
the satellites listed in the legend of the figure. The data are from de- 
tectors having several different energy ranges and no spectral cor- 
rections have been made. The general features of the data in these 
energy ranges have been noted previously in the literature. The flux 
increases rapidly with L, goes through a maximum near L = 1.5 and 
then decreases. The decrease is not as rapid as the initial rise and in 
this energy range the flux generally does not decrease monotoni- 
cally 18, 20 for L > 2. Excepting the measurements of Imhof and 
Smith, 32 the flux decreases with increasing energy, indicating a falling 
energy spectrum. 

The clashed line in Fig. 49 is the ratio of the 50-130 MeV proton flux 
measured with Telstar® 1 to the 40-110 MeV proton flux measured with 
Explorer 15. This ratio is a good qualitative index of the energy spec- 
trum near 45 MeV, and in these circumstances the change in this index 
is independent of the absolute calibrations of the instruments. The 
ratio is seen to decrease monotonically as L goes from 1.25 to 1.9, 
indicating, in agreement with the references cited in the previous sub- 
section, a softer spectrum* at higher L. 

12.5 Intensity vs B on L Shells 

In Fig. 50 [parts (a), (b), and (c)] measurements from various 
satellites are compared on the three L shells, 1.3, 1.5, and 1.8. The 



* A softer .spectrum contains a larger fraction of low-energy particles. 
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Fig. 50 — Values of omnidirectional flux on various specified .L-shells, for the 
satellites and energy ranges indicated in the legend, are plotted vs B, in parts 
(a), (b), and (c). Ratios of Telstar® 1 to Explorer 15 measurements are shown 
in part (d). References are given in Table VII. 
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Explorer 15 and Injun 3 measurements have been compared in more 
detail by Valerio. 10 Observe that J decreases monotonically with B 
on all the L shells and the shape of J vs B is very similar for all the 
measurements on the same shell except for the lowest L shell where 
the dependence on the energy response of the detector is most im- 
portant. Information concerning the energy spectrum near 45 MeV 
is contained in the changes in the ratios of the measurements, and in 
these circumstances the changes are independent of the absolute 
calibrations of the instrument. 

To cast more light on the qualitative behavior of the energy spec- 
trum, the ratio of the 50-130 MeV proton flux measured with the 
Telstar® 1 satellite to the 40-110 MeV proton flux measured with Ex- 
plorer 15 has been calculated as a function of B for fixed L. The results 
are plotted in Fig. 50(d). All the ratios increase with increasing B for 
L from 1.2 to 1.9 inclusive. The values of B in the plot cover the range 
from the magnetic equator to a magnetic dipole latitude (A) of about 
30°. The increase in the ratio indicates a spectrum that hardens with in- 
creasing B in the neighborhood of 45 MeV. At L = 1.8 Frcden et al 35 
find a spectrum that hardens with increasing B for proton energies 
between 10 and 35 MeV, but softens with increasing B for proton 
energies above about 55 MeV. Our results suggest that this change in 
behavior cannot have occurred below 50 MeV. 

12.6 The Intensity Near the Top of the Atmsophere 

The position of the 8-protons/cm 2 sec flux contour from the Telstar® 1 
satellite is plotted in B,L coordinates in Fig. 51 (a), together with our 
own extrapolation of the published Injun 3 data 19 to a flux of about 10 
protons/cm 2 sec,* and the 16-proton/cm 2 sec flux contour from Explorer 
4. The purpose of this figure is to test whether or not the altitude 
dependence of contours of constant counting rate at low altitudes is 
consistent with other data. The qualitative agreement of the results 
plotted in Fig. 51(a) is quite good, especially for L < 1.8, where the 
atmosphere is controlling. A number of effects may contribute to the 
divergence of the results for L > 1.8. Among them are: temporal ef- 
fects, this region of the belt is shown to be subject to temporal varia- 
tions in Section X; instrumental effects, the instruments arc near their 
threshold sensitivities in a region of magnetic space in which the en- 
ergy spectrum may be anomalous; and biases in the fitting procedure, 



* Valerio 19 states that his fit* (and therefore his Fig. 8) are not intended to 
represent the data accurately at low altitude. 
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examination of residuals give some indication of a slight bias in the 
fitted function in this region. 

It is difficult to get direct insight into the altitude dependence 
from a B,L plot, so the values of B have been transformed into mini- 
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Fig. 51 — Comparison of isoflux contours obtained from three satellites near 
the top of the atmosphere. Part (a) B, L coordinates, part (b) minimum altitude 
(near the South American magnetic anomoly). References are given in Table VII. 
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mum altitude by using the procedure mentioned in Section XI. The 
minimum altitudes are plotted against L in Fig. 51 (b). It is character- 
istic of all three bodies of data that the minimum in the minimum 
altitude curve does not occur at minimum L. 

It is tempting to consider whether the lower altitude of the Explorer 
4 points, coupled with the lower low-energy threshold and high flux asso- 
ciated with the Explorer 4 measurements, might imply that the exo- 
sphere was less dense when the Explorer 4 measurements were made. 
However, the uncertainty in the position of the Telstai® contour (see 
Section XI) is so large that the use of this figure to refute the hy- 
pothesis that the atmosphere contracted 23 between 1958 (^ solar 
maximum), when the Explorer 4 measurements were made, and 1962, 
when the Telstai-® data were taken is precluded, even if one were pre- 
pared to overlook the possibility that the energy spectrum at these low 
altitudes is anomalous 30 and consequently that the calculated geometric 
factors of the instruments may be in substantial error near the cutoff. 

12.7 Equatorial Pitch Angle Distribution 

The solid curves in Fig. 52(a) represent the equatorial pitch angle 
distributions, at various values of L, calculated from (8) and the co- 
efficients in Table V. When these are compared with the equatorial 
pitch angle distributions obtained from the Injun 3 data, 10 which have 
been replottcd as the dashed curves in Fig. 52(a), they are found to 
be very similar in shape, although the Telstar® curves are a trifle 
flatter. This would be anticipated from the previous discussion of the 
tendency of the energy spectrum of protons with energies near 45 MeV 
to harden at high values of B. The shape of the distributions are, how- 
ever, appreciably different from those derived by Lenchek and Singer 37 
from consideration of possible injection and loss mechanisms. This 
may be seen in Fig. 52(b) which contains the present results as the 
solid lines, and the results of Lenchek and Singer 37 as the dashed lines. 

12.8 Other Bodies of Data 

A portion of the considerable body of relevent high-energy proton 
data, some of which does not meet the requirements for inclusion in 
Table VII, is noted here. The earliest measurements of proton intensi- 
ties were made on Explorers 1 and 3 by Van Allen. 38 His historic esti- 
mate of ^ 2 X 10 4 protons/cm 2 sec with energies >40 MeV at the heart 
of the inner belt (x = 0, L ^ 1.56) has been substantiated by all the 
measurements reported to date. In particular, the high-energy proton 



1416 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 19(57 



800 



700 - 



•v 500 - 



400 - 



9 

2 


1 II 


[MOl 


1 


1 




i 




(a) 


i 


- 






«.._ 
















^Z^f^^^^O 


? 








2 




o^ 








\\ 






L^ 


VJVI.70 


? 

5 

2 


1 1 1 1 


L = 1.20 

I I ! 


I 

j 1.25 

[pa 

i i 


1.30 

1 


1.40 

1 






.80 
\l.90 

\ ~ 
12.00 - 




0.9 1.0 



Fig. 52 — Unidirectional flux vs x on various L-shells. The solid curves are de- 
rived from the HTB coefficients using (8). The dashed curves in part (a) are 
Injun 3 results (from Valerio 1 '-' Fig. 8). The dashed curves in part (b) are the 
results of the theoretical calculations of Lenehek and Singer" taken from their 
Fig. 10 and arbitrarily normalized to reasonable values of j at x = 0. 
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measurements made in the inner belt by Explorers 6, 12, and 14 and 
Pioneers 3 and 4, have been noted by Frank et al 3!) to agree with each 
other and with those on Explorers 1 and 3. Reference to some measure- 
ments made with ballistic probes may be found in the article by Freden 
et al. 3S 

XIII. QUO VADIS 

The mathematical model which has evolved along the lines summa- 
rized in Section IV has provided a very satisfactory representation of 
the high-energy proton data from the Tehtar® 1 satellite, as discussed 
in both statistical and physical terms in Sections VI through XII. It is 
appropriate to consider how this work might be extended. 

13.1 Further Improvements within the Present Scheme 

The final fit of Model II has a mean square error which is less than 
twice the variance to be expected on the assumption of a Poisson dis- 
tribution of the count data. Some of this excess is surely due to "ex- 
perimental error." However, one might seek some additional improve- 
ment by the addition of more parameters to the fitting function as 
indicated in Model III of Section 4.7. Such fits, carried out on an 
approximately 1000- point selected data set, will almost surely lead to 
a reduction in the mean square residuals because of the increased free- 
dom the additional parameters provide. However, as noted in Section 
4.7, preliminary work with Model III has not led to a really substan- 
tial improvement, either statistically or aesthetically as judged by plots 
of the residuals. 

Additionally, one might try to improve further on the representative- 
ness of the sample by simple iteration. Using the HTB fit to Model II 
to determine new x,L cells, another sample might be selected and fitted. 
The very small differences between the Model-I CB fit and the Model-I 
(or II) HTB fit do not suggest that this would be fruitful in the present 
case. If the preliminary fit used for determining the x boundaries of 
the cells were a poorer fit, iteration would clearly be worthwhile. 

A further extension of the procedure for designating representative 
cells would involve the development of a two-dimensional version of 
the basic idea and procedure outlined in Section 7.1. Specifically, one 
would try to define approximately 1000 x,L cells within each of which 
the preliminary fit to y(x, L) has the same range. In the present case, 
the anticipated gain from this refinement did not seem to justify the 
practical difficulties. However, a practical, well-defined algorithm for 
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such a process in several dimensions simultaneously might prove very 
useful. 

13.2 Another Approach to the Model 

All the models presented so far are of the form 

y(x,L) = A(L)-b(x;e,(L)), (18) 

where A(L) represents the variation in intensity along the magnetic 
equator and b(x; e f (L)) represents the variation with x on an L-shcll. 
The {e,(L)} adjust the nature of the dependence on x, as a function 
of L. This approach arises from the L-shell orientation of the adiabatic 
theory of trapped particle motion. 

Alternatively, one might focus attention on the shape of y as a func- 
tion of L at constant x, rather than on y as a function of x at constant 
L. It is shown in Fig. 19(a) and discussed in Section 7.2 that y(x, L) 
as a function of L for fixed x forms a simple nesting set of curves at 
successive values of x. This is a consequence of the monotonic decrease 
of y with x at any fixed L. With this orientation, a model might be 
expressed as: 

y(x,L) = F(L;p,(x)), (19) 

where F{L) is the shape of a constant-re section, whose parameters, 
the {pi}, are expressed as functions of x. Although this approach would 
not contain the L-shell orientation of the particle motion explicitly, it 
seems to offer veiy significant practical possibilities. 

13.3 Full Data Utilization 

In the two-dimensional fits that were carried out, only a selected set 
of data were used, either chosen at random within a set of narrow, 
contiguous L-slices, as in the fit of Section VI, or chosen on the basis 
of a preliminary fit to the data as in Section VII. All the data were 
examined by residual plots and mean square residual measures of the 
fits, but only a small part of the data were actually used in determin- 
ing the values of the fitting parameters. With this procedure, informa- 
tion is clearly being lost that could be used to "better" determine the 
function. 

Several methods have been applied in the past to allow all of an 
existing body of satellite data to influence the mathematical descrip- 
tion of that data. The most direct method uses interpolation or smooth- 
ing functions. It is often the case that consecutive satellite observa- 
tions from a particular detector are closely enough spaced to determine 
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the local spatial variation. Under those conditions a sequence of data 
points can be averaged or fitted to a local smoothing function. A num- 
ber of points in a sequence may thus be replaced and represented by a 
single point which is determined by them all. The replacement may 
also be made at some particularly convenient coordinate location, for 
example, at one of a fixed set of L or x values on which functional 
fitting may subsequently be carried out. This method has been used 
on the data of Explorer 15, portions of which have been described by 
Mcllwain, 18 Roberts 10 and Brown. 41 

In the context of the high-energy proton data from the Telstar® 1 
satellite, a different but analogous procedure could be used. Rather 
than selecting at random one data point within each of approximately 
1000 x,L cells, all points within a given cell could be used to determine 
a value which would represent the observable at the central point of 
the cell. This might be done by simply averaging the points within the 
cell, but the cell size is large enough so that the x and L dependence 
within the cell generally cannot be neglected. A more representative 
procedure would be to fit the points within an x,L cell with a local 
smoothing function. This function can be the same function with which 
the finally selected data values would be fitted across the complete 
range of x,L space (see Appendix B.7). Although in the present case 
the average number of points per cell is about 40, in many cells the 
number of points is fewer than the number of coefficients of the Model 
II function, and some coefficient constraint would be required. This is 
not a substantial objection, however, since the function is only being 
used for smoothing and does not need to be capable of elaborate varia- 
tion over an x,L cell. 

A procedure of this kind greatly reduces the chance that members 
of a final 1000-point set will be nonreprescntative and acknowledges 
the experimental weight of adjacent observations in fixing the values 
of the set. Accordingly, one would expect a reduction in the mean 
square residuals overall the data, from a fit to such a smoothed sample. 

The procedure of smoothing within a cell could be used with larger 
x, L cells (with more points per cell) to define a point set smaller than 
1000. It can of course also be used with much larger bodies of data 
up to a maximum of 1000 points per cell with the existing computer 
program. 

13.4 Extension to Other Cases 

There are very evident values in being able to communicate the 
essence of a large body of data in terms of a mathematical model with 
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a small number of coefficients. This is very effectively accomplished by 
the present empirical representation of the Telstar® 1 high-energy 
protons, but the model is very specialized. As previously noted, includ- 
ing a wider range of space such as that explored by the Telstar® 2 
satellite requires modification of the function. Characterizing the pro- 
ton distribution for substantially lower energy protons may well re- 
quire functions outside the generality of even Model III. Treating 
electrons in almost any region of space requires treating time as well 
as position variables because a complete set of measurements of the 
spatial distribution of the particles cannot readily be obtained in a 
time short compared with significant time variations. 

No single formulation yet exists which is capable of coping in a use- 
ful way with the range of measurements of particles trapped in the 
magnetic field of the earth. However, the success of the present formu- 
lation as it has been evolved and the general methods that have been 
developed gives us confidence that other and more complicated cases 
can be treated. 

XIV. SUMMARY AND CONCLUSIONS 

This section provides a summary, with references, for the entire 
document including the appendices. 

14.1 General A ccomplishment 

The main accomplishment is the development of a relatively simple 
(empirical) mathematical model which gives a statistically accurate 
representation of the spatial distribution of high-energy protons meas- 
ured with the Telstai"® 1 satellite. 

14.2 The Data 

14.2.1 Space and Time Coverage (Sections I and II) 

The data were acquired between July 1902 and February 19G3 within 
the region of space bounded by 1.09 R e ^ R g 1.95 R. and ^ X ^ 58°. 
Inside these boundaries good temporal and spatial coverage were 
achieved. 

14.2.2 Energy Range and Instrumental Sensitivity (Appendix A) 

The nominal energy interval of the detector is 50 < E < 130 MeV 
and its nominal geometric factor is 0.1431S:Sa cm2 ster - Tne m " 
strument is effectively omnidirectional and the lower threshold of 
sensitivity is «1 proton/cm 2 sec. 
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14.2.3 Telemetry (Section II) 

Each observation consisted of the number of counts registered in 11 
seconds. With this was associated the time at which the telemetry was 
received, and auxilliary information. 

14.3 The Models 

14.3.1 Coordinate System (Section III) 

Each model relates the omnidirectional intensity of high-energy pro- 
tons to a two-dimensional magnetic space whose coordinates, x,L, de- 
rive from a mapping of the earth's main magnetic field onto an axially 
symmetric dipole field through the adiabatic invariants of the particle 
motion. 

14.3.2 General Form and Properties (Section IV) 

The models have the form of a product, A(L) -G(x,L), in which 
the first term expresses the equatorial intensity as a function of L, and 
the second term describes the diminishment of intensity, as a function 
of increasing .r, for fixed L. The functional expressions for G (exclud- 
ing G'") transform in closed form to equivalent pitch angle distribu- 
tions. 

14.3.3 Specializations (Sections IV and IX) 

Retrospectively, all the models may be considered to be specializa- 
tions of Model III, but historically the two-dimensional models evolved 
from a series of one-dimensional fits on L-slices. These fits led to the 
L-slice model which was then generalized empirically to the two-di- 
mensional Model I. Model I was in turn specialized to Model II to 
overcome some statistical fnonlincarities and high correlations) and 
interpretive difficulties encountered with Model I. 

14.4 Fitting 

14.4.1 Criterion (Section III and Appendix B) 

The least squares criterion was used in deriving estimates of the 8 
(or 9 or 10) coefficients required by the models to fit the data. 

14.4.2 Scale (Section III and Appendix B) 

To stabilize the variance of the observations, the models have been 
fitted to the square root of the observed counting rate. 
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14.4.3 Sampling (Sections 6.1, 6.9, 7.1 and Appendix B.3) 
Coefficients of Models I and II were estimated by fitting samples 
containing about 1000 of the nearly 80,000 available observations. 
Sampling is necessary to avoid exaggerating the importance of regions 
of x,L space where data are abundant, and also for compatability with 
existing computer programs. A method of sample selection based on a 
preliminary fit has been developed to provide a good overall represen- 
tation of the data. Before selecting the sample, the data were parti- 
tioned to remove instrumental effects and outliers identified by study- 
ing residuals from preliminary fits. 

14.5 Quality oj Fit 

14.5.1 Criteria oj Judgment (Sections VI to IX and Appendices B 
andC) 

Judgments regarding the quality of fit were largely based on graph- 
ical studies of residuals, the behavior of the fit at the boundaries of 
the radiation belt and various statistical measures. Residuals (equal 
to observed minus fitted), on the square root scale, were particularly 
useful as sensitive indicators of the quality and nature of the fit. 

14.5.2 Comparisons Among Models (Sections V and IX) 

The L-slice fits give good one-dimensional representations of very 
limited regions of data. Both the standard errors of the coefficients and 
the correlations among coefficients are high compared to the corre- 
sponding measures derived from the two-dimensional fits. The fits of 
Models I and II to the 960-point HTB sample are practically equiva- 
lent. However, Model II is superior in the following respects: one less 
coefficient is required, standard errors are uniformly smaller, correla- 
tions among the coefficients are uniformly smaller, the index of non- 
linearity is veiy much smaller, and more of its coefficients have a phys- 
ical meaning. 

14.5.3 Coordinates (Sections VI and VII) 

Plots of residuals vs x, L, time, etc. indicate the general adequacy 
of x,L coordinates for the organization of the data. 

14.5.4 Quantitative Measures (Sections VII, VIII, IX, and Appendices 
B and C) 

Typically, the fits account for nearly 99 percent of the variability 
about the data mean. The mean square error of fit is about 1£ times 
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as large as would be anticipated on the basis of assumed Poisson sta- 
tistics. Even in the worst of quite small spatial regions, the mean 
square residual does not exceed 2\ times the Poisson-based prediction. 
Probability plotting procedures indicate that the residuals are closely 
normally distributed and lead to an estimate of the variance which is 
about twice the Poisson-based prediction. 

14.5.5 General Limitations (Appendix C) 

Statistical examination of all the data, categorized in x,L cells de- 
fined from a preliminary fit, indicates that it is unlikely that the fit 
given by the present model could be significantly improved by any 
simple modification based on x,L coordinates alone. 

14.6 Numerical Values of Fitted Coefficients, Standard Errors, etc. 

14.6.1 L-Slices (Section V) 

Coefficient values and other statistics for four L-slices appear in 
Table II, and values of coefficients for a large number of L-slices are 
shown in Figs. 8 to 10. 

14.6.2 Models I and II (Sections VI to IX, also Sections V, XI, and 
XII) 

Model II is the preferred model. Coefficients, standard errors, cor- 
relations, and other summary analysis-of-variance statistics appear in 
Table IV for Model I and Table V for Model II. The coefficient func- 
tions: (?) square root of average counting rate, y(x,L); {ii) square 
root of average equatorial counting rate, A(L); and (Hi) position of 
cutoff, x c (L); are well-determined and applicable values, standard 
errors, and correlations appear in Figs. 19 and 24 for y(x, L) (and 
Figs. 45 to 47 for the flux) ; Figs. 8, 11. 21. 22, and 30 for A(L) ; and 
Figs. 9, 12. 22, 23, and 30 for x,.(L) (and Fig. 44 for altitude). 

M .7 Some Physical Results 

M.7.1 Flux Maps (Section XII) 

Flux maps are given in B,L and R,\ coordinates and as J,B contours 
for constant L, based on the fitted model and using a calibration of the 
detector assuming certain single-component energy spectra. Neglecting 
uncertainties of calibration, the relative fluxes have a standard error of 
about 2 percent. The value of the maximum flux is (S.Tl^g) X 
10 3 protons/cm 2 sec at L = 1.46 on the magnetic equator. 
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14.7.2 The Cutoff (Section XI) 

The minimum geographical altitude corresponding to the fitted cut- 
off function was computed. This altitude varies as a function of L and 
has a value of about 270 km at the magnetic equator at L = L<> = 1.13 
and a minimum of about 160 km at L = 1.6. The shape of this L de- 
pendence suggests that, the interaction between the protons and the 
residual atmosphere is of major importance in determining the cutoff 
for values of L less than 1.9. For larger L values, the loss mechanism 
determining the cutoff is of different origin. 

14.7.3 Temporal Effects (Section X) 

The general spatial distribution of high-energy protons is very stable 
in time over the period covered by the present data; however, using 
residuals as a sensitive indicator, we find two temporal effects that are 
distinguishable from instrumental effects. Firstly, there appears to be 
an increase in the flux in the 1.9 < L < 2.2 region during the 30-day 
period starting about day 280, 1962. This increase varied from about 
5 to 90 percent depending on both x and L. Secondly, there is an indica- 
tion of a qualitative increase in the altitude of the cutoff over the pe- 
riod of the observations. The present results indicate that any diurnal 
variability of the earth's magnetic field would have an upper limit of 
0.003 Gauss at L ^ 1.5. 

14.7.4 Comparison with Other Experiments and Theory (Section XII) 
The absolute fluxes measured in this experiment agree well (within 

a factor of two) with other extensive experimental measurements, but 
the present values are in general slightly lower. Spatial distribution of 
the flux agrees very well with other measurements but differs appreci- 
ably from published theoretical calculations. 

14.8 Extensions (Sections XIII, IV, and Appendix B) 

The methods developed in this work have lead to a very satisfactory 
representation of the high-energy proton data from the Telstai® 1 
satellite. 

With the better methods of utilizing data and selecting samples 
noted in this paper, and with more general functional forms (some 
approaches to which have been indicated), it should be possible to rep- 
resent the radiation intensity for other more extensive and less "well- 
behaved" bodies of data than the one treated here. Most aspects of 
the statistical methods developed are generally applicable to problems 
of modeling data mathematically. 
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appendix A. 

The Instrument 

Energetic electrons and protons were measured on the Telstar® 1 
satellite by a group of detectors in all of which the sensitive element 
was a phosphorous-diffused silicon diode specially developed for such 
particle measurements. 7 The active volume of the device is the disk- 
shaped space-charge region of the diode under reverse bias. For the 
detector measuring protons with energies above 50 MeV, the reverse 
bias was approximately 100 volts, the space-charge region was approx- 
imately 2.8 mm in diameter and 0.39 mm thick, and the diode was 
shielded by about 12 mm of aluminum over a solid angle of 2tt and 
somewhat more than 12 mm of aluminum equivalent over the remain- 
ing hemisphere (see Fig. 53). 

The thickness of the space-charge region of the detector was meas- 
ured with protons from a cyclotron. A calculation of the path-length 
distribution for unscattered particles in the space-charge region and in 
the surrounding shielding materials has been made. These calculated 
results have been combined with range-energy information, and the 
properties of the associated electronic circuits, to give the geometric 
factor of the instrument, g(E), as a function of the energy, E, of pro- 
tons incident on the spacecraft. The geometric factor varies with the 
reverse bias voltage and the temperature of the detector, both of which 
affect the effective thickness of the active volume of the diode. Fig. 54 
is a graph of g(E) vs E for a bias voltage of -97.5 volts and a temp- 
erature of 20°C, the nominal operating conditions of the instrument. 
Note that protons with energies below 50 MeV were not detected. 

The geometry of the detector and shield is only approximately omni- 
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Fig. 53 — The instrument. 

directional. However, the satellite was spin stabilized, the symmetry 
axis of the detector was nearly perpendicular to the spin axis of the 
satellite, and the telemetered counting rate was an average over at 
least 15 revolutions of the satellite. This averaging process tends to 
remove any directionality inherent in the detector geometry. A sensi- 
tive analysis noted in Section 7.10 failed to show any directional de- 
pendence in the data. 
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Fig. 54 — Dependence of geometric factor on energy of protons incident on the 
shielding. 
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For a differential energy spectrum N(E), where N(E)dE is the number 
of protons with energies between E and E + dE, the average geometric 
factor, g(E u E 2 ), of the detector for particles with energies between 
E { and E-, is defined by 



S(E lf E 2 ) = 



f g(E)N(E) dE 
f ' N(E) dE 

Jr. 



(20) 



The function </(50 MeV, Z? 2 ) has been evaluated numerically for various 
values of E 2 and forms of N(E). The values of 0(50 MeV, 130 MeV) 
are plotted in Fig. 55 as a function of n for the single-component power- 
law spectrum N(E) <x E~ n , and also as a function of E for the single- 
component exponential spectrum N(E) « exp(—E/E ). It may be 
seen from the figure that 0(50 MeV, 130 MeV) varies by less than 
G percent from 0.143 cm 2 ster for < n < 7.5 and 10 MeV < E < 
90 MeV. These ranges of n and E include most experimentally de- 
termined values by a comfortable margin. 3, °' 29,34,35 The omnidirectional 
flux, J(E U E 2 ), of protons with energies between E, and E 2 is given by 



J(E U E 2 ) -- 



4ttF 2 



g(E 1} E 2 



(21) 



where Y is the counting rate of the detector. In the body of this paper, 
the values E x = 50 MeV, E 2 = 130 MeV and 

g = 0(50 MeV, 130 MeV) = 0.143 cm 2 ster (22) 

are used. The flux J(50 MeV, 130 MeV) is designated simply by J, 
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and the counting rate to flux conversion is considered to be independent 
of the proton energy spectrum. 

While the relative value of g shows a variation of less than 6 percent 
for the wide range of single-component energy spectra noted above, 
the absolute value of g is less well specified. Variations in the ambient 
temperature and reverse bias voltage may change the effective geo- 
metric factor by as much as 25 percent. The difficulty of dealing with 
the complexities in shielding geometry, caused by embedding the 
instrument in the spacecraft, introduces additional uncertainties in 
the absolute value of g. These uncertainties are in the range of —25 
to +50 percent. 

No provision was made for recalibrating the detector once the 
satellite was in orbit. However, the evidence, which is discussed in 
Section X, concerning the temporal variations of the proton distribution 
is that neither the detector nor the associated circuit elements were 
substantially affected by the space environment. Instrumental (e.g., 
temperature and bias voltage) effects are often quite different in char- 
acter from temporal changes in the proton belts and may be separated 
from them in many circumstances. It is, of course, possible to postulate 
instrumental effects that will be inextricably confounded with certain 
secular changes that might take place in the proton distribution. 

APPENDIX B. 

Some Statistical Details 
B.l Introduction 

This appendix presents, heuristically, some facts and formulae con- 
cerning the statistical analysis of the data. While a variety of statis- 
tical principles, precepts and procedures were employed as guides, the 
main judgments came from empiricism, scientific intuition and com- 
mon sense. Various kinds of plots of residuals, used informally, have 
been of key importance, both for evaluation and for suggestion. 

Simply stated, the objective was to produce a statistically accurate 
analytical description of the intensity distribution of high-energy pro- 
tons in space surrounding the earth. The process of analysis involved 
the empirical evolution of a mathematical model, in interaction with 
the application of fitting and evaluative techniques. The data source 
and processing have been described in Sections II and III. The itera- 
tive and interactive processes of the final stages of model development, 
fitting, data partitioning and data sampling are described in Sections 
IV to IX. 
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Appendix B.2 deals with the basis for use of the square root trans- 
formation, }', of the counting rate data, Y 2 . Appendix B.3 discusses the 
selection of a sub-sample used in fitting. The use of the method of least 
squares in nonlinear model fitting, to estimate unknown coefficients, 
or functions of the coefficients, and their standard errors and correla- 
tions is reviewed briefly in Appendix B.4. Some remarks on construc- 
tion of sums of squares contours, often referred to as confidence re- 
gions, and of indices of local nonlinearity of the model are given in 
Appendix B.5. Appendix B.6 discusses several issues relevant to the 
interpretation of the analysis of variance results. Appendix B.7 de- 
scribes a mode of "smoothing" data within cells, which could have 
been used in conjunction with the sub-sampling procedure. Appendix 
B.8 concerns the technique of probability plotting. 

B.2 The Square Root Transformation 

It appears a reasonable assumption (supported by some empirical 
evidence) that, in the absence of geophysical disturbances, at a fixed 
point in space relative to the earth, the number of counts Z, recorded 
in the detector in 11 seconds, will vary in time according to a Poisson 
distribution, i.e., 

Probability {Z = z) = 6 -~- , z = 0, 1, 2, 3, • • • , (23) 

z\ 

where the parameter of the distribution, v, is the mean value of Z. 

With this statistical model, the average intensity of radiation in the 
region of space measured by the detector is proportional to v, where 
the proportionality factor depends on the counter geometry and effi- 
ciency. The objective is to develop a function which describes how v 
varies in space, based on observations of the quantity Z at different 
positions in the satellite orbit. 

For the Poisson distribution, the variance of Z is also v, i.e., the 
average of the squared deviations, (Z — v) 2 , is v. Thus, as the value of 
v changes, the variance of the associated random variable Z also 
changes. Hence, the scatter of Z about its average value will be differ- 
ent in different regions of space as the average intensity fluctuates. 

Working with the experimental data on the scale of Z has two draw- 
backs. Firstly, if one fitted a mathematical model to the data using a 
least squares criterion, the different observations would have variable 
weight, which would require appropriate, troublesome, allowance in 
the fitting procedure. Secondly, graphical judgment of the adequacy of 
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any particular fit would be difficult because, of the variable scatter of 
the data about a fitted function in different regions. 

Thus, the square root transformation, Y, of the counting rate 
(Y 2 = Z/ll counts per sec) was used to "stabilize" the variance and 
the model-fitting procedure employed unweighted nonlinear least 
squares on Y (but with some data weighting as discussed in Section 
7.1 and Appendix B.3). 

Heuristically, consider the linear Taylor's expansion of Z about v 

VZ± V~v + { ^ 1 --- . (24) 

Then, the variance of y/~Z is approximately 

Var (VZ) = (-T/=) 2 Var (Z - v) + • • • . (25) 

If 

Var (Z - v) « v, (26) 

then 

Var(VZ) cc j;- j, (27) 

that is, Var (V-Z) would be approximately a constant. 

Discussions of this transformation are given by Bartlett 42 and 
Anscombe. 43 If the distribution is in fact Poisson, then Anscombe shows 
that the average value of \/~Z is approximately 

V; - 8 ^ " .28,' ' 
while the variance of y/Z is, asymptotically, 

1 ii . 3 4- XL . 
I V + & + 32, 2 + 

Again for the Poisson distribution, Bartlett gives exact values of the 
dependence of the variance of VZ on v, summarized in the following: 

p:\0\ 0.5 J 1 ! 2 1 3 1 4 I 6 1 9 1 15 
Var VZ : |0 10 .3 10 10 .402 10 .390 10 .340 10 .306 10 .276 10 .263 10 .256 

For a Poisson distribution, a transfo rmation of the form s/Z +1/2 
or VZ + 3/8 or (VZ + VZ + 1 - 1) will improve the variance 
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stabilization at smaller v values. In the present application, such a 
modification would have appeared physically artificial and incon- 
venient. Moreover, the actual variance of the observations exceeds 
the Poisson variance (see Appendix C) and the "correction" was thus 
felt to be unwarranted. Some response to the (empirically defined) 
variance instability remaining after the square root transformation 
was made in the form of some weighting in the data selection (described 
in Section 7.1 and Appendix B.3). 

Of course, if one wished to adopt the assumption of a Poisson dis- 
tribution as an absolute basis for procedure, instead of as a guide, then 
one might choose to use maximum likelihood to estimate the coeffi- 
cients of the model. This would mean developing a procedure for de- 
termining values of the coefficients [of the function v(x,L)] which 
would maximize 

II e" ( " L) Wx, L)] ! /z\. 

observations 

111 the present case, a general program for nonlinear least squares was 
available while a procedure for Poisson likelihood maximization would 
need to be evolved. Apart from this practical consideration, however, 
it seemed more robust to use the Poisson assumption as a guide to 
developing an appropriate transformation preliminary to fitting by 
least squares. The point is that the square root trans formation will 
effect an approximate variance stabilization not only when the variance 
is equal to the mean (as in the case of the Poisson distribution) but 
also when, more generally, the variance is proportional to the mean. 
Empirical vindication of this caution is given in Appendix C. More- 
over, the least squares approach enables the approximate, statistical 
interpretation of results using familiar procedures from linear multiple 
regression methods. 

The present analysis is based on the quantity }', where Y 2 = count- 
ing rate = Z/ll counts per see. Thus, if in fact Z were a Poisson vari- 
able, 

Var(V)^ (£)(!) = 0.023, (28) 

as a reasonable approximation. When the average counting rate ex- 
ceeds 1/11, this value of 0.023 is a lower bound on the variance of Y, 
even with the Poisson assumption. Moreover, there are many other 
possible sources of intrinsic variability and experimental error in this 
situation. 
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A further benefit which one might expect from the square root trans- 
formation in this circumstance is that the distribution of residuals 
would tend to be more symmetric and more nearly normal (Gaussian) . 

Some empirical properties of this square root transformation in the 
present body of data are given in Appendix C. 

B.3 Sample Selection 

As a practical requirement, the available multivariable, multicoeffi- 
cient, nonlinear least squares fitting program could operate with a 
maximum of 1000 data points. Hence, the 41,135 HTB observations 
needed to be sampled or condensed at a 1 in 40 ratio. 

As in all real sampling or experimental design situations, many com- 
peting criteria and practical difficulties were relevant. Perhaps the 
overriding point, explicitly understood here (and probably true in most 
actual model fitting problems) , is that the model which was being de- 
veloped was not the "truth" but was really just a smoothing function 
which one wanted to fit well over a wide region of space. Thus, it was 
not appropriate to think of estimating the model coefficients, say, so 
as to optimize their apparent (indicated) statistical reliability, nor 
would it be appropriate to try to use all the available data in an 
equally weighted manner, since accidents of orbital position and in- 
strumental behavior would have too great an effect on the distribution 
of data points. 

The procedure developed for the present use is outlined in Section 
7.1, with pertinent remarks also in Section 13.3 and Appendix B.7. 

The method of Section 7.1 yielded 960 observations to which the 
model was then fitted using unweighted least squares. The 960 sampled 
observations were selected so as to be roughly speaking, "widely 
spaced," the metric being change in average counting rate. Thus, the 
challenge of fitting the 960-point sample, as measured by sum of 
squares of residuals, is greater, on a per-observation basis, then would 
be that of fitting the entire body of 41,135 HTB observations, very 
many of which are quite close together. The "model bias" difficulties of 
the entire body of data are concentrated in the sample. The statistical 
fluctuation would be approximately the same, on a per observation 
basis, in the sample as in the whole body of data. 

B.4 Estimation Procedure 

The unspecified coefficients of the models defined in Section IV were 
estimated so as to minimize the sum of squares of deviations between 
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the observed }' and fitted y, for the sample array of data. The itera- 
tive, multivariable, multicoefficient, nonlinear least squares fitting 
was executed using a computer program due to Huyett and Wilk," 
based on a procedure outlined by Wilk 43 (see also Lundbcrg, Wilk and 
Huyett I. '"• l: 

The classical statistical properties of least squares estimation, 
namely unbiased estimates with minimum variance, apply in the case 
of statistically uncorrected observations having equal variances and 
with the coefficients to be estimated occurring linearly in the model 
(see, for example, Wilks"). In the present case, even with the square 
root transformation, the observations do not have equal variances but, 
for practical purposes, the weighting implied by the selection proced- 
ure (see .Section 7.1) compensates adequately. The model is, however, 
quite nonlinear in the coefficients. Still, one hopes that the attractive 
statistical properties of linear least squares cany over approximately 
to the nonlinear case because, in small enough neighborhoods, non- 
linear functions can be linearly approximated. (An index for measur- 
ing model non linearity is described in Appendix B.5.) In any case, 
the least squares criterion is geometrically appealing and primitively 
meaningful. 

Among the by-products of the fitting procedure, applied to the par- 
ticular array of data in .v,L space, are approximate values for the 
standard errors of the estimated coefficients, a matrix of approximate 
pairwise correlation coefficients for the estimated coefficients, an anal- 
ysis-of-variance table giving the sum of squares accounted for and 
not accounted for by the fitted model, a list of residuals (equal to 
observed minus fitted), and various plots. 

The least squares estimates of single- valued functions of the coeffi- 
cients, such as A(L), x e (L), or y(x, L) are simply the same functions 
of the estimates of the coefficients (since least squares is an invariant 
process). Approximate variances and correlations of functions of the 
coefficients may be derived as follows: If 6' = (0, , • • • , P ) denotes 
the coefficients of the model, and 6 their estimates, then the approximate 
covariance of the estimates g(d) and h(6) of the functions g(6) and 
h(6) is 

Covariance (g(6), h(6) 

= Cov (g(8), h(6)) 

= Statistical average of \(g(0) - g(6))(h(0) - h(6))\ 



1434 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1907 



~ y y f dg(g) 

_ i i L 9d { ; J 



.^1 Gov (MO 

9 _ do,- Je 



-stWlPffl^ 



+ 2 £ fa^l [Wl tVar (10 Var («0]V„ , (29) 

where p„ is the correlation of 0, and d, . The formula for the approx- 
imate variance of g(d) is then just a specialization of the above, putting 

9 = h - 

Some associated facts and issues are worth mentioning here. First, 

the approximate statistical correlations p,,- of the estimated coeffi- 
cients of the model, or of functions of these, depend on (i) the distribu- 
tion of the sample in x,L space, (ii) the values of the coefficients and 
(fit) the nature of the mathematical model ; but do not depend on the 
actual adequacy or appropriateness of the fit. Similarly, the approximate 
standard errors of estimates are each made up as a product of which 
one term depends upon the square root of the mean square of the 
residuals of fit and the other depends only on the same factors as do 
the pa . Second, the various statistical measures, such as standard 
errors of estimated coefficients which are obtained from the fit to the 
960-point HTB sample are, in a narrow statistical sense, conservative 
because they refer to the sample only and do not make allowance for 
the fact that the fitted model does indeed fit very well to the entire 
body of 41,135 HTB data. Thus, if statistical fluctuations were the 
only factor in the uncertainty of the estimates, one might further 
reduce t his uncerta inty by some factor, roughly approximated by 
6 « V41, 135/960. This view of statistical uncertainty does not, 
however, give appropriate weight to the "model bias", which will not 
be eliminated by any number of observations. Third, all the summary 
statistical measures, which are referred to as standard errors, correla- 
tions, confidence regions, etc., should be used and interpreted in a 
data analytic way, i.e., as indicating facets of the body of data and the 
adequacy of its description by the model and analysis — rather than in 
terms of some supposedly "true" model or hypothesis which one is 
trying to evaluate in probabilistic terms. 

B.5 Sums of Squares Contours, "Confidence Regions" and 
Nonlinearity Indices 

The models of Section IV are defined up to the values of the un- 
specified coefficients. Any set of values for these coefficients may be 
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said to provide a "fit" to the 960-point sample of data. Thence one can 
define a sum of squares function of the set of coefficients as 

SS (coefficients) = £ (observed - "fitted") 2 , (30) 

which will take on various (positive) values as one varies the values 
of the coefficients. In the space of the coefficients there exist then, in 
principle, contours of this "sum of squares" function. 

While standard errors provide information on reasonable allowances 
for the estimate of a single parameter in the light of the fit of the 
model to the actual body of data, they do not carry any information 
on the joint statistical properties of the estimates. A reasonable (ro- 
bust and primitive) indication of joint statistical behavior is provided 
by these "sum of squares" contours in coefficient space. 

In the case of models in which the unknown coefficients occur lin- 
early, these contours are a family of ellipsoids defined by certain sim- 
ple quadratic functions of the coefficients. The orientation and shape 
of this family of ellipsoids indicate the interdependence of the esti- 
mates of the coefficients in the light of the data, and show which 
coefficients are well-determined and which poorly. However, the in- 
terpretive value depends heavily on geometrical appreciation and, for 
more than a few coefficients, high-dimensional representation cannot 
be achieved directly. 

The ellipsoid (even in the linear case) is not defined, in general, by 
its one-dimensional projections. (The standard error of a coefficient 
estimate is half the length of the projection of the unit ellipsoid of 
the family onto the coefficient axis.) But, as a matter of simple geo- 
metrical fact, all pairs of two-dimensional projections do uniquely de- 
fine the ellipsoid. Thus, one practical means of a complete graphical 
representation of the high-dimensional ellipsoid is in terms of all possi- 
ble pairwise planar projections. 

For the case of linear models, on the basis of a series of assumptions 
— namely that the differences between the model and the observations 
are due to statistical fluctuations which are normally and independ- 
ently distributed all with zero mean and the same variance — some may 
choose the abstract probabilistic interpretation of these ellipsoids as 
"confidence regions" (see, for example, Wilks 16 ). If this interpretation 
is used, it is necessaiy that the distinctions and relationships between 
the joint, pairwise and marginal confidence coefficients and regions or 
intervals be understood. Details will not be provided here. Briefly, if 
a nine-dimensional ellipsoid were specified to have a confidence coeffi- 
cient of /? 9 , then any two-dimensional projection would have a con- 
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fidence coefficient of 2} interpreted marginally. The relation between 
ft, and p-, is indicated by the following: 



0.13 0.00 

0.25 0.95 

0.50 0.984 

0.75 0.997 

0.90 0.9994 

0.95 0.99995 
In the present case, the model is nonlinear and the fluctuations are 
not normal. Contours of the sums of squares function as a function of 
the coefficients can, in principle, be obtained for a given body of data 
and will not be ellipsoids. In practice, however, obtaining these con- 
tours is so laborious as to be virtually impossible. 

However, one may consider a linear (planar) approximation to the 
nonlinear model in the neighborhood of the least squares estimates of 
the coefficients and thence obtain expressions for a family of ellipsoids 
which may be reasonably good approximations to contours of the sums 
of squares function. An index of the effective nonlinearity of the 
model is the nonconstancy of the sums of squares of residuals on these 
ellipsoids and this can be normalized by division by the value of the 
minimum sum of squares. Such measures are presented and discussed 
in Sections VIII and IX. 

Given that the linear approximation is adequate, the nonnormality 
of the observations should not deter those who seek (and who believe 
in) the general probabilistic confidence interpretation since the statis- 
tical process is likely very robust. 

Sections VIII and IX contain specific examples of some of the pair- 
wise projections of these "approximations to sum of squares contours." 
Specifically, the size of the 9-dimcnsional ellipsoid was such that, if 
all the statistical assumptions applied, a joint 0.99 confidence coeffi- 
cient could be attached. Since a complete set of pairwise projections 
for nine coefficients involves 36 ellipses only a few are shown. As a 
summary indicator of the nature and behavior of these ellipses the 
quantity 

a = (sign of p)-(l - Vl - P a ) (31) 
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is tabulated (in Tables IV and V), where p is the correlation of the pair 
of coefficients involved. The value of 1 — \<x\ is the ratio of the area of 
the actual ellipse to that of the largest ellipse which could be inscribed 
in the rectangle formed by the horizontal and vertical tangents to the 
actual ellipse (see Wilk 48 ). The range of a is— l^a^l and large 
values of |a| (say above 0.75) corresponds to narrow ellipses with major 
axis oblique to the coordinate axes, and represent situations of high 
interdependence of the coefficient estimates. 

B.6 The Analysis of Variance 

The analysis of variance provides a summary description of the 
apportionment of the "variability" of a body of data in the light of 
the model employed for analysis, where variability is defined in terms 
of sum of squares. 

Given n observations, one may visualize an /^-dimensional observation 
space, whose coordinates represent the possible values of each of the 
n observations. The data are then represented by a fixed point in this 
space. 

The model, having p unspecified coefficients, implies certain functional 
relationships amongst the coordinates of the observation space. Thus 
the model effectively defines a constraining "surface" of p dimensions, 
and each point on this surface corresponds to some set of values of 
the unspecified coefficients of the model. The least squares estimate 
of the coefficients corresponds to that point on the constraining surface 
which is closest to the actual data point. If the coefficients in the model 
occur linearly then the constraining surface is a hyperplane which 
ordinarily, by definition of the observations, contains the origin, and, 
if the model includes a constant term, also contains the equiangular 
line (corresponding to the mean). 

The squared distance of the data point to the origin is then the total 
sum of squares, 2-i ^? < while its shortest squared distance to the 
constraining surface is the error or residual sum of squares, associated 
with lack of fit. The difference between these may be termed the model 
sum of squares and, for linear models, this is actually the squared 
distance from the least squares estimates point to the origin.* If a 
constant term is included in a linear model, then the model sum of 
squares may be further decomposed additively in terms of the squared 



* In the linear case, the model sum of squares is easily computed directly as 
the squared length of the projection onto the hyperplane of the line joining the 
data point, and the origin. This fact is used in the present iterative computer pro- 
gram in checking convergence. 



1438 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1907 

perpendicular distance (call it D\) of the least squares estimate point 
to the equiangular line and the squared distance (call it D\) along the 
equiangular line, from the foot of that perpendicular, to the origin. 
This latter quantity D\ is usually termed the sum of squares due to 
the mean. The squared distance of the above-defined point on the 
equiangular line to the data point is called the corrected total sum 
of squares, £ (Y t - Y) 2 and is just 2 Y\ - 1% . The ratio of the 
squared length D\ to the corrected total sum of squares is denned as 
the squared multiple correlation, R 2 , and often used as a measure 
of accomplishment of a model. It is easy to show that R 2 denned above 
is equal to 

sum of squares for error 

total corrected sum of squares 
This latter quantity is computable even when the model is nonlinear 
and/or docs not contain a constant term. 

One may define contours of sums of squares of residuals in the con- 
straining surface as the loci of the intersections with the surface of 
given radii from the observation point. In the event that the constrain- 
ing surface was a hyperplane, which would be true if the unspecified 
coefficients in the model occur linearly, then these loci (or contours) 
would be a family of p-dimensional spheres. For nonlinear models, 
this will be approximately true for a sufficiently small neighborhood 
of the least squares point. 

The particular form of the model, in regard to the unspecified co- 
efficients, defines a coordinate system within the constraining surface. 
Three cases are worth distinguishing. First, the constraining surface 
is a hyperplane and the coefficients are linear. Second, the surface is 
a hyperplane but its coordinates are nonlinear. The second case may 
be reduced to the first by appropriately transforming the coefficient 
coordinate system. Third, the surface is nonlinear. In this case one 
can approximate the surface by a hyperplane in a small neighborhood. 
Thus, in a sufficiently small neighborhood, the situation can be re- 
garded as linear. 

The approximately or exactly linear coordinates implied by the 
model will in general be nonorthogonal. Thus, the representation of 
the spherical (exact or approximate) contours in an orthogonal co- 
ordinate system for the coefficients yields a family of ellipsoids. In the 
sense of measuring lack of fit by sums of squares between fitted and 
observed values, these contours in coefficient space constitute sets 
whose members are "equidistant" from the data point. 
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B.7 A Procedure for Smoothing in Cells 

In Sections 7.1, 13.3 and Appendix B.3, discussion of why and how 
to sample and possibly "smooth" the data has been given. One specific 
practical possibility is now described. 

Suppose one has a preliminary fit of the model, represented by 
g(Wi ; 0), where 9' = {9 X , 9, , ■ ■ ■ , 9 V ) are the fitted coefficient values 
and Wi denotes the independent variables. Suppose this preliminary 
fit is used to partition the space of the independent variables (here 
x and L) into some approximation of equirange cells, as described 
earlier. As argued in Section XIII, it may be profitable to "smooth" 
the data in each cell so as to yield a value generally representative 
of all the observations in that cell, instead of using a random selection 
from the cell. 

A sensible smoothing function for each cell is, clearly, the model 
g(w; 9). A simple procedure is, for each cell separately, to carry out 
one stage of linear adjustment, doing the linear least squares regression 
of { F, — g(w, ; 9) \ on dg/dd { \§ , • • • , dg/dd p \§ , to obtain the regression 
coefficients 8' = (8 X , • • • , 8 P ), for that particular cell. Then the smooth- 
ing function for that cell would be g(w; 9") where 9 = 9 + 8. A rep- 
resentative "smoothed observation" for that cell might then be the 
quantity g(w; 9), where w is, say, the mid-point of the cell. 

This process permits each cell, overall, to determine a single value 
to represent it in the entire fitting process and diminishes the chance 
that a random selection from a cell may be unnecessarily nonrepresenta- 
tive of that cell behavior. 

If one had wished to fit to all the available data, then the smoothed 
cell values would be weighted in proportion to the number of data 
in the cell. In the present case, this was deliberately not done. 

The goodness of fit of a model to smoothed cell values, not dif- 
ferentially weighted, cannot be statistically judged directly from the 
analysis of variance since the residuals are no longer individually 
statistically comparable and the mean square residual is not an estimate 
of the error variance of the observations. However, the fitted model 
can be assessed by functions of its residuals from the original data 
(or a sample thereof). 

B.8 Probability Plotting 

The techniques of probability plotting are useful for data analysis 
in a wide variety of circumstances. (See Wilk and Gnanadesikan 17 for 
a general discussion of probability plotting techniques.) For instance, 
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in the present work, plots of residuals against various variables have 
provided invaluable guidance, but one is also interested in the dis- 
tributional behavior per se of the collection of residuals. As presented 
in Section 8.1, normal and half-normal probability plots have been 
used for this purpose. 

The rationale for such probability plots is roughly as follows: If one 
draws a random sample of size n from a population which is normally 
distributed with mean /* and variance o- 2 then the ordered observations 
would be expected to approximate, roughly, to a linear function, /x + 
<rZi{n), of appropriate "representative" values zdn) from a standard 
normal (/x = 0, a 2 = 1) distribution. Thence a plot of the ordered ob- 
servations against the Zi(n) would tend to be linear, with intercept 
approximately /x and slope approximately u. For the representative 
value, Zi{n), corresponding to the ith ordered observation, one can use 
the standard normal quantile for the proportion (i— h)/n. 

This plotting technique displays the. individual observations in a 
sample graphically and does so against a backdrop such that the ex- 
istence of outliers and asymmetry, as well as other distributional prop- 
erties, are sensitively indicated. Of course such plots are usually profit- 
ably supplemented by others that order or partition the data according 
to information extraneous to the responses themselves. 

We expect the mean of the residuals, Y - y, in the present study (see 
Section 8.1) to lie near 0. Also we expect that their variances will be 
approximately the same, since that is the purpose of the square root 
transformation. As a further benefit of the square root transformation 
we expect that the distribution of the residuals will tend to be sym- 
metric and to approach normality; thence the present application of 
normal probability plotting of the residuals. The fact that these resid- 
uals are not entirely statistically independent — since they derive from 
a commonly estimated fitted function — is a minor issue since the num- 
ber of observations is so much larger than the number of fitted coeffi- 
cients. 

Half-normal probability plotting employs the ordered absolute re- 
siduals plotted against standard half-normal (standard normal folded 
at 0) distribution quantiles. Such a plot eliminates any symmetry-type 
information but provides an incisive focus in bringing together on the 
plot the largest departures from fit. 

Probability plots can provide very sensitive indications of distribu- 
tional peculiarities especially in regard to "overly" large values. Some- 
times the indications are of little practical interest, such as minor 
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lumps which one can sec in Fig. 29. but in other regards, such as in 
estimating an "intrinsic" error standard deviation, the plots may per- 
mit a good judgment on how to discount or correct for apparently 
aberrant values which might otherwise have an undue influence, say, 
on mean square error. 

Error standard deviations may be estimated from normal or half- 
normal probability plots as the "slope of the linear configuration." 
Typically, it will not be relevant to make a great show of objectivity 
in this process since the declared purpose is to permit an informal dis- 
counting of unexpected distributional peculiarities. Thus, in Fig. 29, 
one takes the slope as defined essentially by the bulk of residuals, 
ignoring the few largest. 

appendix c 

Statistical Measures Over All the HTB Data 

This appendix presents various statistical measures over all the 
41,135 HTB data. These measures concern the fit of Models I and II 
and the partition of the x,L space (as described in Section 7.1 and 
Appendix B.3) into 1034 cells of which 813 were nonempty of obser- 
vations. The partition is such that the range of y within cells is rela- 
tively small. For each cell, two functions are used: (i) The mean 
square deviation (MSI) I defined as 

MSD (u) = -1— f) (u, - uf, (32) 

ii i j = i 

where the cell has n observations and u, denotes some function of a 
cell observation, e.g., F, or Y 2 , and ii is the mean of the ?/, in the cell; 
(ii) The mean square residual (MSR) defined as 

MSR(I') =\ Eir,-!/,) 2 , (33) 

n , = i 

where (/, is the fitted value (from Model I or III corresponding to the 
observed Yj. 

C.l Empirical Justification of Square Root Transformation 

Figs. 56 and 57 show plots of MSD(T 2 ) versus the cell mean of Y 2 
and MSD (Y) versus the cell mean of Y, respectively. It is seen that 
MSD(F 2 ) shows a distinct and major dependence on the average 
value of the counting rate, Y-, while MSD (7) does not show syste- 
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raatic increase relative to the average value of Y, except, as expected, 
in the close neighborhood of zero counting rate. 

A more detailed analsis of the results of Fig. 56 indicates that the 
dependence of MSD(Y 2 ) on cell mean of Y 2 is somewhat curvilinear 
having larger slope for larger Y 2 values. This curvilinearity is very 
likely mainly due to the mode of definition of the x,L cells. The 
procedure used tends to produce cells which are "too large" in regions 
where the counting rate is also large, thus leading to an apparent 
extra increase in MSD(Y 2 ) with Y 2 . At all values, however, the 
dependence of MSD(Y 2 ) on Y 2 is greater than the slope 0.09 ( = 1/11) 
which would be associated with the Poisson distribution. The em- 
pirically observed slope varies from about 0.15, based on small values, 
to 0.3, based on large values of the MSD(Y 2 ). 

These results suggest that one cannot hope to achieve, by means of 
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Fig. 56 — Cell MSD (F 2 ) vs cell mean of Y s for the x,L cells denned in Sec- 
tion 7.1 and Appendix B. 
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Fig. 57 — Cell MSD (F) vs cell mean of Y for the x,L cells defined in Section 
7.1 and Appendix B. 

any fitted model based on x,L coordinates, on the scale of Y, a mean 
square residual (error) as low as 0.023 which is associated with the 
Poisson assumption. 

Although the Poisson assumption provided a useful stimulus toward 
a profitable transformation of the data, these results confirm that 
it would have been unwise to have tied oneself too closely to the 
assumption as a complete basis for analysis, as for instance in basing 
the fit on maximization of the Poisson likelihood function (see 
Appendix B.2). Possible sources of variability and error in the data, 
beyond Poisson fluctuations in counts, have been discussed elsewhere 
in this paper. 

C.2 Determination of Weights 

The sample selection procedure involved "weighting" the 813 non- 
empty cells by selecting 2, 3/2, or 1 observation per cell. The observed 
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MSD (F) were classified into three groups defined by: ^ MSD ^ 
0.013; 0.013 < MSD ^ 0.02; 0.02 < MSD. The x,L coordinates of the 
midpoints of cells so identified are shown in Fig. 58. The actual assign- 
ment of weights was based on applying contiguity considerations to this 
plot. 

C.3 Analysis of Variance Over All the HTB Data 

Table VIII summarizes the analysis of variance over all the 41,135 
HTB data. The table covers the fit of Models I and II to all the data, 
using the estimated coefficients (see Tables IV and V) from the fit 
to the 960-point sample. Also, one can regard the collection of 
averages of the Y values in each of the 813 nonempty cells as pro- 
viding a fit depending on 813 fitted quantities. The corresponding 
"error" (cell deviations) is the pooled cell MSD (Y). Finally, the 
residuals of the fit of Model I (or II) can be "fitted" by 813 cell 
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Fig. 58 — Positions of centers of regions in x,L space having certain ranges of 
cell MSD (Y). The ranges are indicated in the legend. 
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Table VIII — Analysis of Variance over all the HTB Data 
(41,135 Points Minus 226 Outliers). 



Due to 


d.f. 


Sum of squares 


Mean square 


Total (41,135-226) 


40,909 


230,267.45 




Mean 


1 


115,755.39 




Corrected total 


40,908 


114,512.07 




Model I residuals 


40,900 


1,411.3 


0.0345 


Model II residuals 


40,901 


1,419.6 


0.0347 


Cell deviations 


40,096 


1,541.4 


0.03S4 


Cell dev. of Model I res. 


40,085 


1,167.0 


0.0291 


Celldev. of Model II res. 


40,086 


1,166.9 


0.0291 



Multiple R- value 



Model I 
Model II 



0.98S 
0.9SS 



averages of the residuals, leading to an "error" which is the pooled 
cell MSD(y - y), i.e., due to the cell deviations of the Model I (or 
II) residuals. 

The fits to all the data provided by Models I and II are equally 
good, as was true for the 960-point sample. The mean square resid- 
uals over all the data (about 0.035) is lower than the value (about 
0.036) obtained for the sample even though the fit of the model was 
determined by the sample. This bears out the expectation (see Ap- 
pendix B.3) that the mode of selection of the sample is such that the 
sample was harder to fit on a per-observation basis -than the entire 
body of data. 

The cell means provide overall a poorer actual fit than Model I or II, 
and allowing for the number of fitted coefficients, the mean square for 
cell deviations exceeds that for the models by about 12 percent. 

Fitting cell means to the model residuals yields an additional sub- 
stantial reduction in the sum of squares of the model residuals and a 
mean square of about 0.029, which is some 17 percent lower than the 
value for the models. If in fact the models gave an "unbiased" fit every- 
where, then one would expect that the values of pooled MSR(Y) and 
pooled MSD(l r — y) would be nearly the same. The excess of the 
former is due mainly to systematic inadequacies of the fit (see 
Appendix C.4). 

The value 0.029 represents virtually a lower bound on the achiev- 
able 'mean square error' for this body of data. Despite its downward 
bias from the substantial number of 'zero counting rate' observations, 
it exceeds the 'Poisson' variance of 0.023 by about 25 percent. This 
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excess is probably due to a combination of factors, including incom- 
plete elimination of temperature and bias voltage instrumental effects, 
as discussed further below. 

The 'improvement' of the MSD( Y - y) over MSR cannot be taken 
to mean that some smooth "simple" adjustment of the model based on 
x,L coordinates might be found so as to yield similar improvement. 
Some of the bias apparently associated with x,L coordinates in dif- 
ferent regions may be due to artifactual association with temporal, 
instrumental, or other small effects and such corrections could not be 
made overall in terms of a "simple" x,L dependence. 

C.4 Analysis of the Excess Variation 

A study of plots of cell MSD(Y) against each of y, x, and L 
indicates that large MSD values occur mainly in the 1.2 < L < 1.7 
region, at high average counting rates. This excess is due largely to 
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Fig. 59 — Cell MSR (F) vs central value of y' for cell. 
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Fig. 60 — Positions of centers of regions in x, L space having certain ranges of 
cell MSRO r )/MSD(r). The ranges are indicated in the legend. (Plotted 
points are mid-points for the cells. Points appearing to the right of the boundary 
R = 2.0 R, represent cells which have data only in one corner.) 

the hybrid mode of x,L cell formation, in which the L-slices were 
equal length intervals, while within each //-slice, the x segments were 
chosen to have equal increments of y. Thus, at L values where y is 
large, the x,L cells will tend to have a larger y range. 

The tendency of MSD to rise with cell average counting rate is 
not mirrored by cell MSR behaviour. As shown by Fig. 59, the level 
of MSR is not dependent on ?/ except, as expected, for those cells 
where the counting rate is near zero. Roughly speaking, the average 
level of cell MSR for y values away from zero is about 0.04, in agree- 
ment with the probability plot estimate of Section 8.1. Of course, Fig. 
59 shows both smaller ordinate values and less dependence on the 
abscissa values than the comparable plot of Fig. 57. 

The relation of cell MSR to cell MSD is partially indicated in Fig. 
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60, showing positions in x,L space of various ranges of the value of 
the ratio MSR/MSD. One sees that MSR tends to exceed MSD along 
the "outside" of the data region. The excess along the R = 2 R, 
boundary is due mainly to model bias or inadequacy. The excess at 
high L-high x is probably associated with temporal effects. The large 
ratios along the left edge of the data, which is the cutoff region, is 
likely a reflection of deficiency of the function. The excess of MSR 
over MSD is associated in the main with small y values. 

Fig. 61 shows cell mean square deviations of residuals, MSD( Y —if) , 
plotted against y. This plot shows less scatter (most noticably for 
MSD(F - y) > 0.075) than that of Fig. 59, and a lower average level 
of MSD ( Y - y ) for y > 0, as expected. The high values of MSD (Y — y) 
are not related to y as such but rather, as other plots show, with 
"extra fluctuations" in the 1.2 < L < 1.7 region. This is probably 
associated with the coarse HTB data partition which does not com- 
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pletely take care of the temperature and bias voltage instrumental 
effects. 
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